Prejudgment control method and system for main shaft fault of numerical control machining center

Abstract

The invention relates to the technical field of numerical control machine tool predictive maintenance and health management, and discloses a numerical control machining center spindle fault pre-judgment control method and system, and the method comprises the steps: obtaining current data and temperature data in the start-stop process of a spindle motor, recognizing the type of a data segment based on a start-stop instruction timestamp, and generating a labeling data set; performing Hilbert transformation on the three-phase current data of each stage and extracting a current envelope to generate a frequency band sub-sequence set; calculating a mahalanobis distance based on the health domain boundary, and generating an impact resistance degradation index; and predicting the growth trend of the topological sensitivity gradient by using a damage accumulation criterion and a crack propagation formula, calculating the residual start-stop times reaching a critical threshold, and generating a residual life prediction result. Quantitative modeling and residual life prediction of the shock resistance degradation process of the main shaft are achieved, and sudden failure of the main shaft under the condition of frequent start and stop is avoided.

Description

Technical Field

[0001] This invention relates to the field of predictive maintenance and health management technology for CNC machine tools, and more specifically, to a method and system for predicting and controlling spindle failures in CNC machining centers. Background Technology

[0002] The high-speed CNC machining center spindle in the flexible manufacturing unit uses ceramic hybrid bearings and can reach speeds of up to 30,000 rpm. It is frequently started and stopped according to workpiece switching needs, with more than 50 start-stop cycles per shift. The CNC machining center spindle is driven by a three-phase asynchronous motor. Temperature sensors are installed in key parts of the spindle to monitor the spindle temperature, and current sensors are installed at the motor input to collect three-phase current signals.

[0003] Under frequent start-stop conditions, the multi-band topology of the spindle motor current envelope undergoes drastic changes during each start-stop transient process. Simultaneously, rapid fluctuations in temperature and current create transient thermo-electric shocks. A healthy spindle exhibits rapid recovery of its current envelope topology and a stable response to transient shocks. However, a degraded spindle experiences slow recovery of its high-frequency topology, permanent deformation, and significantly increased sensitivity to transient shocks, ultimately leading to a loss of shock resistance and sudden failure.

[0004] Existing technologies rely on independent analysis of current topology characteristics or transient impact intensity for fault prediction. However, the degradation of spindle shock resistance under frequent start-stop conditions is a result of the combined effects of multi-band topology stability and the cumulative effects of transient thermo-electric shocks. Existing technologies have the following problems: single-band topology analysis cannot capture degradation signs where weak high-frequency signals are masked by low-frequency components; steady-state topology monitoring ignores the dynamic response characteristics of the topology during transient processes and does not quantify the sensitivity changes of the topology to cumulative impacts, leading to inaccurate predictions of spindle fatigue resistance degradation under start-stop impacts and sudden system failure under cumulative damage. Summary of the Invention

[0005] This invention provides a method and system for predicting and controlling spindle faults in CNC machining centers, solving the technical problems in related technologies where independent analysis of temperature or current cannot comprehensively reflect impact intensity; independent analysis of topological features or impact intensity cannot reflect the coupling effect between the two; low-frequency components mask weak high-frequency signals in full-band topology analysis; and steady-state topology monitoring ignores the dynamic response of transient processes.

[0006] This invention provides a method for predicting and controlling spindle faults in CNC machining centers, comprising:

[0007] Acquire current and temperature data during the start-up and shutdown process of the spindle motor, identify the data segment type based on the start-up and shutdown command timestamps, and generate a labeled dataset;

[0008] The three-phase current data at each stage are subjected to Hilbert transform and the current envelope is extracted. The data is then decomposed into multiple frequency band subsequences using a multi-resolution filter bank to generate a set of frequency band subsequences.

[0009] Construct a topological space for each frequency band sub-sequence and calculate the persistent homology features to generate a multi-frequency band topological feature matrix;

[0010] Based on the multi-band topological feature matrix, high-frequency anomaly indicators and topological sensitivity gradient tensors are calculated to generate topological dynamic response features.

[0011] Based on the temperature rise rate and the coefficient of variation of current amplitude, a two-dimensional thermo-electric phase space trajectory is constructed, the geometric characteristics of the trajectory are calculated, and transient impact characteristics are generated.

[0012] The transient impact intensity is quantified using the rainflow counting algorithm. A joint feature space of the topological sensitivity gradient tensor and the equivalent impact stress amplitude is constructed, and a joint feature space trajectory is generated.

[0013] Calculate the curvature of the joint feature space trajectory and the frequency of abrupt events of high-frequency anomaly indicators, calculate the Mahalanobis distance based on the health domain boundary, and generate an index of shock resistance degradation.

[0014] The growth trend of topology sensitivity gradient is predicted by using the damage accumulation criterion and crack propagation formula, the remaining number of start-stop cycles before reaching the critical threshold is calculated, and the remaining lifetime prediction results are generated.

[0015] Furthermore, the steps of constructing the topological space for each frequency band sub-sequence and calculating the persistent cohomology features include:

[0016] For high-frequency segment subsequences, Takens' delayed embedding theorem is used to construct a delayed embedding phase space, which reconstructs the one-dimensional time series into a set of points in an m-dimensional phase space;

[0017] The optimal delay time is determined by the mutual information method, and the delay value corresponding to the first local minimum value of the mutual information curve is selected.

[0018] The optimal embedding dimension is determined by the pseudo-nearest neighbor method. The embedding dimension corresponding to the pseudo-nearest neighbor ratio is the optimal value when it drops below 5%.

[0019] Topological connections at different scales are made to the point set in phase space, and simple complex sequences are formed by the Vietoris-Rips complex construction method;

[0020] Calculate the generation and disappearance of topological features at various scales to generate persistent homology barcodes;

[0021] The statistical features of the barcode length are extracted as topological invariants, including mean, variance, maximum value, and entropy, which constitute the topological feature vector.

[0022] Furthermore, in the step of calculating the high-frequency anomaly index based on the multi-band topological feature matrix, the calculation formula for the high-frequency anomaly index is the energy ratio of the high-frequency topological feature vector relative to the energy ratio of the mid-frequency topological feature vector and the low-frequency topological feature vector, that is, the square of the Euclidean norm of the high-frequency topological feature vector divided by the sum of the squares of the Euclidean norms of the mid-frequency and low-frequency topological feature vectors.

[0023] Furthermore, the topology sensitivity gradient tensor is the difference gradient matrix between the topological feature vectors of each frequency band before and after the transient, which contains the difference between the topological feature vectors of the three frequency bands: low frequency, mid frequency, and high frequency. The norm of the gradient tensor reflects the overall magnitude of the topological structure change, and the condition number reflects the degree of non-uniformity of the topological change between different frequency bands.

[0024] Furthermore, the steps for constructing the thermo-electric two-dimensional phase space trajectory based on the temperature rise rate and the coefficient of variation of current amplitude include:

[0025] The first derivative of the temperature time series is calculated using the central difference method as the rate of temperature rise.

[0026] The ratio of the standard deviation to the mean of the current amplitude time series is calculated as the coefficient of variation;

[0027] The temperature rise rate and the coefficient of variation of current were respectively normalized using Z-score.

[0028] Two-dimensional phase space trajectories are constructed based on standardized thermal and electrical state variables;

[0029] The maximum curvature and trajectory length of the phase space trajectory are calculated as geometric features.

[0030] Furthermore, the method for calculating the maximum curvature of the trajectory is as follows: calculate the first and second derivatives of the thermal state variables and electrical state variables of the phase space trajectory respectively, obtain the curvature value at each moment according to the curvature calculation formula, and take the maximum value as the maximum curvature of the trajectory; when the trajectory velocity is less than the preset tolerance, the curvature at that moment is set to zero.

[0031] Furthermore, the steps for quantifying transient impact intensity using a rainflow counting algorithm include:

[0032] The curvature time series of the thermo-electric coupling phase space trajectory is used as the equivalent stress time history;

[0033] The peaks and valleys of the curvature time series are extracted to form an extreme value series;

[0034] Identify complete cycles and calculate the curvature amplitude for each cycle;

[0035] Count the number of cycles for different amplitude ranges;

[0036] The equivalent impact stress amplitude for a single start-stop cycle is calculated by weighting the amplitude of each cycle by combining the trajectory length and transient duration.

[0037] Furthermore, the steps for calculating the Mahalanobis distance based on the health domain boundary include:

[0038] Collect joint feature data from at least 100 start-stop processes under healthy spindle conditions;

[0039] Calculate the mean vector and covariance matrix of the health data;

[0040] Robust covariance estimation is performed using the minimum covariance determinant method.

[0041] The boundary of the health domain is defined as a hyperellipsoid whose Mahalanobis distance is equal to the critical value of the chi-square distribution;

[0042] Calculate the Mahalanobis distance of the current joint eigenvectors relative to the health domain mean and covariance matrix.

[0043] Furthermore, the steps for predicting the growth trend of the topology sensitivity gradient using the damage accumulation criterion and crack propagation formula, and calculating the remaining number of start-stop cycles before reaching the critical threshold, include:

[0044] The cumulative fatigue damage degree of the historical start-up and shutdown process was calculated based on the Palmgren-Miner linear damage accumulation criterion.

[0045] Predicting the growth trend of the norm of the topological sensitivity gradient tensor with the number of start-stop cycles using the Paris crack propagation formula;

[0046] The remaining number of start-stop cycles to reach the critical threshold for the topological sensitivity gradient is calculated using numerical integration.

[0047] The average fatigue life at the current stress level and the remaining number of start-stop cycles based on cumulative damage are calculated based on the SN curve.

[0048] The smaller value between the topology sensitivity prediction and the cumulative damage prediction is taken as the final prediction result of the remaining start-stop count.

[0049] This invention provides a spindle fault prediction and control system for CNC machining centers, comprising:

[0050] The data acquisition module is used to acquire current and temperature data during the start-up and shutdown process of the spindle motor;

[0051] The multi-band decomposition module is used to perform multi-resolution filtering and decomposition on the current envelope;

[0052] The topology feature extraction module is used to calculate the continuous coherence topology features of each frequency band;

[0053] The dynamic response analysis module is used to calculate the topological sensitivity gradient tensor and high-frequency anomaly indicators;

[0054] Thermoelectric coupling analysis module is used to construct thermo-electric two-dimensional phase space trajectories and quantify impact intensity;

[0055] The joint feature analysis module is used to construct a joint feature space and calculate Mahalanobis distance;

[0056] The lifetime prediction module is used to predict the remaining number of start-stop cycles based on damage accumulation and crack propagation.

[0057] The decision output module is used to generate early warning signals and optimization strategies.

[0058] The beneficial effects of this invention are as follows:

[0059] This invention decomposes the current envelope into three frequency band subsequences—low-frequency, mid-frequency, and high-frequency—using a multi-resolution filter bank. It independently calculates the continuous homology topological features for each frequency band, overcoming the factor that low-frequency components mask weak high-frequency signals in full-band topology analysis. This allows it to capture subtle topological changes caused by high-frequency degradation symptoms such as microcracks in ceramic bearings.

[0060] This invention generates a topology sensitivity gradient tensor by calculating the differential gradient of the topology feature vectors of each frequency band before and after the transient, thereby quantifying the dynamic response characteristics of the spindle system topology during the transient process. This overcomes the factor that steady-state topology monitoring ignores the dynamic response during the transient process, and thus can identify the degradation of the topology recovery capability of the spindle under start-stop transient impact.

[0061] This invention achieves quantitative assessment of transient thermo-electric impact intensity by constructing a two-dimensional phase space trajectory of thermo-electricity and using a rainflow counting algorithm to quantify the equivalent impact stress amplitude of each start-up and shutdown, thus overcoming the limitations of independently analyzing temperature or current to comprehensively reflect impact intensity.

[0062] This invention constructs a joint feature space of topological sensitivity gradient tensor and equivalent impact stress amplitude, and calculates the Mahalanobis distance from the joint space trajectory to the boundary of the healthy domain. This enables the coupled analysis of the dynamic response characteristics of the topological structure and the transient impact intensity, overcoming the inability of independent analysis of topological features or impact intensity to reflect the coupling effect between the two. Thus, it can accurately assess the sensitivity changes of the spindle system to impact.

[0063] This invention utilizes the Palmgren-Miner linear damage accumulation criterion to quantify the cumulative fatigue damage during historical start-up and shutdown processes, and combines it with the Paris crack propagation formula to predict the growth trend of the topological sensitivity gradient. This enables quantitative modeling and remaining life prediction of the spindle's impact resistance degradation process, thereby providing accurate early warning before the spindle loses its impact resistance under the combined effects of multi-band topological degradation and cumulative fatigue damage, and avoiding sudden failure of the spindle under frequent start-up and shutdown conditions. Attached Figure Description

[0064] Figure 1 This is a flowchart of a spindle fault prediction and control method for CNC machining centers according to the present invention;

[0065] Figure 2 This is a dynamic piecewise linear diagram of the multi-frequency band decomposition of the transient current envelope during startup, as presented in this invention.

[0066] Figure 3 This is a bar chart comparing the energy of multi-band topological features before and after the transient state of this invention.

[0067] Figure 4 This is a scatter plot of the thermal-electric coupling phase space trajectory and transient impact evolution of the present invention;

[0068] Figure 5 This is an area diagram of the cumulative fatigue damage evolution curve of the spindle according to the present invention;

[0069] Figure 6 This is a hybrid Y-axis plot of the topology sensitivity gradient and remaining lifetime prediction dual-axis trend of the present invention. Detailed Implementation

[0070] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, some features described in the examples may be combined in other examples.

[0071] At least one embodiment of the present invention discloses a method for predicting and controlling spindle faults in a CNC machining center, such as... Figure 1 As shown, it includes the following steps:

[0072] Step 100: Obtain current and temperature data during the start-stop process of the spindle motor, identify the data segment type based on the start-stop command timestamp, and generate a labeled dataset;

[0073] Acquire high-sampling-rate three-phase current data (at least 50kHz) throughout the entire spindle motor start-up and shutdown process, and simultaneously acquire spindle temperature data (at least 10Hz). Utilize the start-up and stop command timestamps provided by the CNC system to segment the acquired data into three phases: a start-up transient phase, a steady-state phase, and a stop transient phase. The start-up transient phase refers to the time from the issuance of the start command to the spindle speed stabilizing; the steady-state phase refers to the time from the issuance of the stop command to the complete stop of the spindle. Generate a labeled dataset containing timestamps, phase labels, and corresponding data.

[0074] It should be noted that the above-mentioned identification of the transient segment at startup means that when the spindle speed reaches more than 95% of the target speed and the speed fluctuation rate is less than 2%, the transient segment is judged to have ended and the steady-state segment has been entered.

[0075] It should be noted that the identification of the above-mentioned transient segment refers to the moment when the spindle speed begins to decrease from the steady-state value as the starting point of the transient segment, and the moment when the speed drops to less than 5% of the target speed as the ending point of the transient segment.

[0076] The spindle of a CNC machining center in a flexible manufacturing unit has accumulated 1200 hours of operation. During the start-up and stop process of machining an aluminum alloy shell workpiece, the CNC system recorded the start-up and stop command timestamp as 2024-11-15 09:23:45, with the target speed set at 18000 rpm. The data acquisition system recorded the three-phase current at a sampling rate of 50kHz and the spindle front bearing temperature at a sampling rate of 10Hz. Based on the speed sensor signal, after a 2.8-second transient start-up phase, the spindle speed reached 17100 rpm (95% of the target speed), and the speed fluctuation rate dropped to 1.7%, indicating the entry into the steady-state phase. The steady-state phase lasted 120 seconds to complete workpiece machining. After the stop command was issued, the spindle speed began to decrease from the steady-state value, dropping to 900 rpm (5% of the target speed) after 3.2 seconds, ending the transient phase.

[0077] Table 1. Data segmentation results for a certain start-stop process:

[0078]

[0079] Step 200: Perform Hilbert transform on the three-phase current data of each stage and extract the current envelope. Decompose it into multiple frequency band subsequences using a multi-resolution filter bank to generate a set of frequency band subsequences.

[0080] Perform a Hilbert transform on the three-phase current time series data obtained in step 100 to extract the current envelope sequence for each stage. The Hilbert transform converts a real signal into a complex analytic signal, and its amplitude is the current envelope. For the extracted current envelope sequence, perform frequency band decomposition using a multi-resolution filter bank, dividing the envelope into three frequency band subsequences: low-frequency, mid-frequency, and high-frequency. The low-frequency band has a cutoff frequency less than 1 kHz, the mid-frequency band ranges from 1 to 5 kHz, and the high-frequency band ranges greater than 5 kHz. Generate a set of frequency band subsequences containing the subsequences for each stage.

[0081] It should be noted that the construction of the above-mentioned multi-resolution filter bank refers to the use of Butterworth filter banks, with low-pass filters for low-frequency bands, band-pass filters for mid-frequency bands, and high-pass filters for high-frequency bands. Each filter is of order 8 to ensure the frequency band separation effect.

[0082] In this embodiment, to improve the accuracy of current envelope extraction, the three-phase current data is preprocessed before performing Hilbert transform to remove power frequency interference and high-frequency noise. Specifically, notch filtering is performed on the three-phase current data with a center frequency of 50Hz and a bandwidth of 2Hz to remove power frequency interference; then low-pass filtering is performed with a cutoff frequency of 20kHz to remove high-frequency noise introduced during sampling.

[0083] After preprocessing the A-phase current data during the startup transient phase, a Hilbert transform was performed to extract the current envelope. The extracted envelope amplitude ranged from 12.5A to 38.7A, with the peak occurring 1.2 seconds after startup. An 8th-order Butterworth filter bank was applied for band decomposition, with the low-frequency cutoff frequency set at 0.8kHz, the mid-frequency bandpass range at 0.8-5.2kHz, and the high-frequency band start frequency at 5.2kHz. After band decomposition, the low-frequency envelope amplitude ranged from 10.2A to 35.4A, the mid-frequency envelope amplitude ranged from 1.8A to 3.2A, and the high-frequency envelope amplitude ranged from 0.5A to 2.1A. The steady-state and shutdown transient phases underwent the same processing procedure.

[0084] Table 2. Frequency band decomposition results of the transient current envelope during startup:

[0085]

[0086] Figure 2 The changes in the current envelope amplitude over time for three frequency bands (low frequency 0-0.8kHz, mid frequency 0.8-5.2kHz, and high frequency >5.2kHz) during the startup transient period (2.8 seconds) are shown. After extracting the envelope using Hilbert transform, frequency band decomposition is performed using an 8th-order Butterworth filter bank.

[0087] Step 300: Construct a topological space for each frequency band sub-sequence and calculate the continuous homology features to generate a multi-frequency band topological feature matrix;

[0088] For the high-frequency segment subsequence generated in step 200, a delayed embedding phase space is constructed, then a sustained homology barcode is calculated, topological invariants are extracted, and a high-frequency topological feature vector is generated. For the mid-frequency and low-frequency segment subsequences, sustained homology barcodes are calculated and topological invariants are extracted, generating mid-frequency and low-frequency topological feature vectors, respectively. The topological feature vectors of each frequency band in the pre-transient, steady-state, and post-transient stages are combined to generate a multi-frequency band topological feature matrix.

[0089] It should be noted that the construction of the aforementioned delayed embedding phase space refers to using Takens' delayed embedding theorem to construct a one-dimensional time series... Reconstructed The set of points in the 3D phase space, reconstructed as the points in the phase space are:

[0090]

[0091] in To delay time, Embedding dimension. Delay time. The embedding dimension is determined using the mutual information method. Determined by the pseudo-nearest neighbor method.

[0092] Furthermore, the aforementioned mutual information method determines the delay time. The specific steps are as follows: Calculate the time series Its delayed sequence Mutual information between For different delay values Starting with one sampling point and gradually increasing it to 1 / 10 of the sequence length, a mutual information curve is plotted. The delay value corresponding to the first local minimum of the mutual information curve is selected as the optimal delay time. The optimal delay time makes and The correlation between them is minimized, ensuring the independence between the dimensions of the embedding space.

[0093] Furthermore, the aforementioned mutual information content The calculation formula is:

[0094] Time series and delayed sequence Discretize the data into several intervals and statistically analyze the joint probability distribution of each interval. and marginal probability distribution , The formula for calculating mutual information is: ,in For discretized interval indexing of time series, The index of the discretized intervals for the delayed sequence is obtained by summing and traversing all interval combinations. The number of discretized intervals is usually the square root of the time series length, rounded down.

[0095] Furthermore, the aforementioned pseudo-nearest neighbor method determines the embedding dimension. The specific steps are as follows: for a given embedding dimension Calculate the nearest neighbor of each point in the phase space, and then increase the embedding dimension to... The process involves checking the distance changes of the original nearest neighbor in high-dimensional space. If the distance increases by more than a preset threshold (usually doubled), it is identified as a false nearest neighbor. The proportion of false nearest neighbors is then calculated. Start by gradually increasing the dimension. When the proportion of false nearest neighbors drops below 5%, the corresponding embedding dimension is the optimal embedding dimension. .

[0096] It should be noted that the calculation of the aforementioned persistent homology barcode refers to performing topological connections on a set of points in phase space at different scales to form a simplex complex sequence, calculating the generation and disappearance of topological features (connected components, holes, cavities) at each scale, and generating a persistent homology barcode. The length of the barcode represents the persistence of the topological features, and the statistical features of the barcode length (mean, variance, maximum value, entropy) are extracted as topological invariants to form a topological feature vector.

[0097] Furthermore, the specific steps for constructing the aforementioned simple complex sequence are as follows: For the point set in the phase space, define a series of increasing distance thresholds. ,in The minimum distance threshold, This is the second distance threshold. The maximum distance threshold, For each distance threshold, the total number of distance thresholds is [number]. ,in For distance threshold sequence indexing, select those with a distance less than [a certain value]. Pairs of points are connected to form edges, creating a simple complex. The distance threshold starts from the minimum pairwise distance of the point set, and gradually increases to the diameter of the point set with a step size of 1 / 100 of the point set diameter. Using the Vietoris-Rips complex construction method, when the pairwise distances of three points are all less than... When a triangle (two-dimensional simplex) is formed, the pairwise distance between the four points is less than 1 / 2. When a tetrahedron (three-dimensional simplex) is formed, a complete sequence of simplex complexes is formed.

[0098] In this embodiment of the application, in order to extract more detailed topological features, when calculating the continuous homology barcode, not only zero-dimensional homology (connected components) and one-dimensional homology (ring structure) are calculated, but also two-dimensional homology (cavity structure) is calculated. The statistical features of the continuous homology barcode in each dimension are extracted respectively, and the features in each dimension are concatenated to form a complete topological feature vector.

[0099] For the high-frequency segment subsequence of the startup transient segment, the mutual information method was applied to determine the optimal delay time as 23 sampling points (corresponding to 0.46 ms), and the pseudo-nearest neighbor method determined the optimal embedding dimension as 5 dimensions. A 5-dimensional delay embedding phase space was constructed, resulting in a point set containing 2351 phase space points. The Vietoris-Rips complex method was used to construct a simple complex sequence, with the distance threshold increasing in steps of 0.015 from the minimum pairwise distance of 0.12 to a point set diameter of 1.98, for a total of 124 scales. Continuous cohomology at each scale was calculated, and 147 zero-dimensional cohomology barcodes, 38 one-dimensional cohomology barcodes, and 12 two-dimensional cohomology barcodes were extracted. The feature values ​​of barcode length in each dimension were statistically analyzed: the mean length of zero-dimensional homogeneous barcodes was 0.186, the variance was 0.043, the maximum value was 0.892, and the entropy was 3.74; the mean length of one-dimensional homogeneous barcodes was 0.067, the variance was 0.021, the maximum value was 0.235, and the entropy was 2.58; the mean length of two-dimensional homogeneous barcodes was 0.032, the variance was 0.008, the maximum value was 0.089, and the entropy was 1.92. These 16 statistical features were concatenated to form a high-frequency topological feature vector. The same process was performed on the mid-frequency and low-frequency subsequences.

[0100] Table 3. Multi-band topology characteristic parameters during startup transient segment:

[0101]

[0102] Step 400: Calculate the high-frequency anomaly index and topological sensitivity gradient tensor based on the multi-band topological feature matrix to generate topological dynamic response features;

[0103] The energy ratio of high-frequency topological eigenvectors to that of mid-frequency and low-frequency topological eigenvectors is calculated to generate high-frequency anomaly indicators for each stage. The differential gradients between the topological eigenvectors of each frequency band before and after the transient event are calculated to generate a multi-band topological sensitivity gradient tensor. The topological sensitivity gradient tensor reflects the dynamic response capability of the principal axis system's topology under transient impact.

[0104] Specifically, the formula for calculating high-frequency anomaly indicators is as follows:

[0105]

[0106] in, , , These are the high-frequency, mid-frequency, and low-frequency topological feature vectors, respectively. The Euclidean norm of a vector is denoted by .

[0107] Furthermore, to ensure the numerical stability of the calculation of high-frequency anomaly indicators, when the denominator... hour( For numerical tolerance, the value is [value to be filled in]. Set the denominator to To avoid division by zero errors, high-frequency abnormal indicators are at this time. It reflects the absolute energy level in the high-frequency band.

[0108] The formula for calculating the topological sensitivity gradient tensor is:

[0109]

[0110] in, , and These are the first and second transients before and after the transient, respectively. The topological feature vector of the frequency band.

[0111] It should be noted that the above calculation of energy proportion refers to using the square of the Euclidean norm of the vector as the energy measure. When the energy proportion of the high-frequency topological feature vector increases significantly, it indicates that there is an abnormal topological change in the high-frequency band, which may correspond to high-frequency degradation symptoms such as microcracks in ceramic bearings.

[0112] In this embodiment, to more accurately quantify the dynamic response characteristics of the topology, the norm and condition number of the topology sensitivity gradient tensor are also calculated. The norm of the topology sensitivity gradient tensor reflects the overall magnitude of the topology change, and the condition number reflects the degree of non-uniformity of the topology change between different frequency bands. The calculation formula is as follows:

[0113]

[0114] in, Denotes the Frobenius norm. and These are the maximum and minimum singular values ​​of the tensor, respectively.

[0115] Furthermore, to ensure the validity of condition number calculation, when the minimum singular value hour( For singular value tolerance, the value is [value to be filled in]. Set the condition number to a predefined upper limit. This indicates that the topology sensitivity gradient tensor exhibits highly ill-conditioned characteristics, meaning that the non-uniformity of topological changes between different frequency bands is extremely severe.

[0116] For this start-up and shutdown process, multi-band topological feature vectors were extracted for three stages: the start-up transient segment (before the transient), the steady-state segment, and the stop transient segment (after the transient). The Euclidean norm of the high-frequency topological feature vector in the start-up transient segment was calculated to be 4.27, the mid-frequency to be 8.94, and the low-frequency to be 15.38, according to the formula... Calculate the topological eigenvectors for each frequency band during the stopping transient segment (after the transient). The high-frequency norm is 5.83, the mid-frequency norm is 9.12, and the low-frequency norm is 15.61. The topological sensitivity gradient before and after the transient was calculated. Among the 16 components of the differential gradient vector in the low-frequency band, the largest component is 0.28; in the mid-frequency band, the largest component is 0.45; and in the high-frequency band, the largest component is 1.68. The topological sensitivity gradient tensor was then constructed. for Matrix, calculate Frobenius norm Singular value decomposition is performed on the tensor to obtain the maximum singular value. and minimum singular value Calculate the condition number .

[0117] Table 4. Energy changes of multi-band topological features before and after the transient:

[0118]

[0119] Figure 3 The changes in the Euclidean norm (energy) of the topological eigenvectors of three frequency bands—the start-up transient band (before the transient) and the stop transient band (after the transient)—were compared and analyzed.

[0120] Step 500: Construct a two-dimensional thermo-electric phase space trajectory based on the temperature rise rate and the coefficient of variation of current amplitude, calculate the geometric characteristics of the trajectory, and generate transient impact characteristics;

[0121] For each transient data segment acquired in step 100, the first derivative of the temperature time series is calculated as the temperature rise rate, which is used as the thermal state variable. The coefficient of variation (the ratio of standard deviation to mean) of the current amplitude time series is calculated as the electrical state variable. Since the thermal state variable has dimensions while the electrical state variable is dimensionless, Z-score standardization is performed on both the thermal and electrical state variables before constructing the two-dimensional phase space to eliminate the influence of dimensional differences. Based on the standardized thermal and electrical state variables, a two-dimensional phase space is constructed, generating a thermal-electrical coupling phase space trajectory. For the phase space trajectory, the maximum curvature and trajectory length are calculated as geometric features to generate a transient impact feature vector.

[0122] It should be noted that the above calculation of the rate of temperature rise refers to the calculation of the temperature time series. The derivative is calculated using the central difference method, and the formula is as follows: ,in This represents the sampling time interval.

[0123] Furthermore, the calculation formula for the Z-score standardization process described above is as follows: For time series... Standardized sequence ,in The mean of the sequence. The standard deviation of the sequence is 0. The standardized sequence has a mean of 0 and a standard deviation of 1, thereby eliminating the difference in dimensions between different variables.

[0124] It should be noted that the calculation of the maximum curvature of the trajectory mentioned above refers to the parameterized representation of the phase space trajectory. (in For thermal state variables, (For electrical state variables), the curvature calculation formula is:

[0125]

[0126] in, , , , These are the first and second derivatives of the corresponding variables, respectively. The maximum curvature of the trajectory is... .

[0127] Furthermore, the first derivative mentioned above is calculated using the central difference method, and the formula is as follows:

[0128] The second derivative is calculated using the three-point difference method, and the formula is as follows: ,in For the sampling time interval, The derivative calculation method is the same, and forward and backward differences are used instead of central difference at the starting and ending points of the trajectory to ensure numerical stability.

[0129] Furthermore, to ensure the numerical stability of curvature calculations, when the denominator... hour( For speed tolerance, the value is [value to be filled in]. The value ) indicates that the trajectory is approximately stationary at that moment. The curvature at that moment is set to 0, meaning that the trajectory is locally approximated as a straight line segment. This situation usually occurs in the stable phase of a transient process.

[0130] It should be noted that the above calculation of trajectory length refers to discretizing the phase space trajectory and calculating the sum of Euclidean distances between adjacent sampling points. The calculation formula is as follows:

[0131]

[0132] in, For sampling point index, The number of trajectory sampling points, For the first Each sampling time.

[0133] For the transient startup phase, temperature data was extracted at a sampling interval of 0.1 Hz (i.e., every 0.1 seconds), totaling 28 sampling points. The temperature rose from an initial value of 42.3℃ to a peak value of 46.7℃. The central difference method was used to calculate the temperature rise rate, with the peak rise rate of 1.85℃ / s occurring 0.9 seconds after startup. The coefficient of variation (COP) of the current amplitude time series was calculated using a sliding window (window length 0.2 seconds) to calculate the local COP at each time point, with a range of 0.08 to 0.34, and the peak occurring 1.3 seconds after startup. Z-score standardization was applied to both state variables: the temperature rise rate was standardized to a mean of 0 and a standard deviation of 1, and the current COP was also standardized. A two-dimensional phase space trajectory containing 28 trajectory points was constructed, with the curvature reaching its maximum at the 18th point (1.8 seconds after startup). / unit distance. Accumulate the Euclidean distances between all adjacent points to obtain the total trajectory length. Unit length.

[0134] Table 5. Characteristics of the thermal-electric phase spatial trajectory during the start-up transient phase:

[0135]

[0136] Step 600: Quantify the transient impact intensity using the rainflow counting algorithm, construct a joint feature space of the topological sensitivity gradient tensor and the equivalent impact stress amplitude, and generate the joint feature space trajectory;

[0137] Based on the maximum curvature and trajectory length calculated in step 500, the equivalent impact stress amplitude for each start-stop cycle is quantified using a rainflow counting algorithm. The topological sensitivity gradient tensor generated in step 400 is combined with the equivalent impact stress amplitude. Since these two types of features have different dimensions and numerical ranges, mean normalization based on the range is performed on each feature before combination, normalizing all features to the range of 0 to 1, thus constructing a joint feature space. In the joint feature space, the joint features of multiple historical start-stop processes are represented as a sequence of trajectory points, generating a joint feature space trajectory.

[0138] The aforementioned rainflow counting algorithm takes as input the curvature time series of the thermo-electric coupled phase space trajectory and outputs as the stress amplitude and mean of each cycle. Specifically, the curvature time series of the thermo-electric coupled phase space trajectory calculated in step 500 is... As an equivalent stress time history, a rainflow counting algorithm is applied to identify each stress cycle, extract the maximum and minimum values ​​of each complete cycle, and calculate the cycle amplitude. For multiple historical start-stop processes, the number of cycles within different amplitude ranges is counted, and the amplitude of each cycle is weighted based on the trajectory length to calculate the equivalent impact stress amplitude of a single start-stop.

[0139] Furthermore, the specific steps for identifying stress cycles using the aforementioned rainflow counting algorithm are as follows: [The text abruptly shifts to a different topic] ...curvature time series... The peaks and troughs are extracted in chronological order to form an extreme value sequence. Starting from the beginning of the sequence, each of the next three adjacent extreme points is examined. The amplitude range formed by the intermediate extreme value and the preceding and following extreme values ​​is calculated. When the amplitude formed by the intermediate extreme value is contained within the amplitudes of the preceding and following extreme values, it is identified as a complete cycle, and the maximum value of the cycle is recorded. and minimum value Calculate the cyclic amplitude The cycle is removed from the extreme value sequence, and the remaining sequence is identified until no new complete cycle can be identified. Finally, all identified cycles are classified and statistically analyzed according to their magnitude.

[0140] The formula for calculating the equivalent impact stress amplitude is:

[0141]

[0142] in, For circular indexes, The number of identified loop types, For the first The number of loop iterations, This represents the total number of loops. For the first Curvature amplitude of a loop For the trajectory length, For reference trajectory length, This is the influence coefficient of trajectory length.

[0143] Furthermore, the length of the aforementioned reference trajectory The trajectory length influence coefficient is determined by the historical average trajectory length of the healthy spindle under standard start-stop conditions. The experimental calibration value is 0.3. The trajectory length influence coefficient reflects the contribution weight of trajectory length to the equivalent impact stress amplitude.

[0144] In this embodiment of the application, to more comprehensively characterize transient impact characteristics, the influence of the duration of the transient process on cyclic fatigue is also considered when calculating the equivalent impact stress amplitude. The modified formula for calculating the equivalent impact stress amplitude is as follows:

[0145]

[0146] in, The actual duration of the transient process. For reference duration, The duration affects the index.

[0147] Furthermore, the aforementioned reference duration The duration impact index was determined by the historical average transient duration of the healthy spindle under standard start-stop conditions. The experimentally calibrated value is 0.5. The duration effect index reflects the nonlinear influence of transient duration on fatigue damage. When the duration of the transient process is extended, the equivalent impact stress amplitude increases according to a power law.

[0148] Figure 4 The evolution of the thermo-electric two-dimensional phase space trajectory within the 2.8-second transient phase of startup is demonstrated.

[0149] Step 700: Calculate the curvature of the joint feature space trajectory and the frequency of abrupt events of high-frequency anomaly indicators, calculate the Mahalanobis distance based on the health domain boundary, and generate an index of shock resistance degradation.

[0150] For the joint feature space trajectory generated in step 600, calculate the local curvature of the trajectory. Statistically analyze the frequency of abrupt events exceeding a preset threshold for high-frequency abnormal indicators calculated in step 400. Based on historical health status data, establish a health domain boundary in the joint feature space. The health domain boundary is defined using the mean and covariance matrix of the health data. Calculate the Mahalanobis distance from the current joint feature space trajectory point to the health domain boundary, generating an index of shock resistance degradation.

[0151] Furthermore, the specific steps for calculating the local curvature of the aforementioned joint feature space trajectory are as follows: The joint feature vectors of the historical multiple start-stop processes are... Arranged in chronological order to form a sequence of trajectory points, among which This is the joint feature vector of the first start-stop process. This is the joint feature vector of the second start-stop process. For the first The joint feature vector of each start-stop process, For the total number of historical start-stop processes, for the th For the nth trajectory point, calculate the local curvature of that point. The specific method is as follows: using the nth trajectory point... Taking one point as the center, and one point before and one point after it to form a three-point sequence. ,in Let be the joint feature vector of the (i-1)th start-stop process. Let be the joint feature vector of the i-th start-stop process. For the joint feature vector of the (i+1)th start-stop process, calculate the vector. and Local curvature is defined as the rate of change of the angle between the two vectors divided by the arc length, and the calculation formula is: ,in For vectors and The included angle, through Calculated.

[0152] Furthermore, the preset thresholds for the aforementioned high-frequency abnormal indicators This was determined through statistical analysis of historical data on the health status axis, specifically the mean of high-frequency abnormal indicators in the health data. Add 3 standard deviations ,Right now The preset threshold follows the 3σ criterion, ensuring that the probability of exceeding the threshold in a healthy state is less than 0.3%, thereby classifying events exceeding the threshold as high-frequency topological anomaly mutation events.

[0153] Furthermore, the specific method for calculating the frequency of the aforementioned mutation events is as follows: within a fixed time window... Internal statistical high-frequency abnormal indicators Exceeding the preset threshold Number of start-stops The formula for calculating the frequency of mutation events is: The time window The duration of the most recent work shift (usually 8 hours) is taken. The mutation event frequency reflects the density of high-frequency topological anomalies per unit time. The higher the frequency, the faster the spindle degrades. When the mutation event frequency exceeds twice the average mutation event frequency under healthy conditions, it is judged to be in the accelerated degradation stage.

[0154] Furthermore, establishing the aforementioned health domain boundary refers to collecting joint feature data from at least 100 start-stop processes of the main axis in a healthy state, and calculating the mean of these feature vectors. Covariance Matrix The boundary of the health domain is defined as the Mahalanobis distance equal to the critical value of the chi-square distribution. The hyperellipsoid, in which For the feature space dimension, With a significance level of 0.01, the health domain boundary contains 99% of the health status data points.

[0155] It should be noted that the calculation of the Mahalanobis distance mentioned above refers to the calculation of the current feature vector in the joint feature space. Relative to the mean of the health domain Covariance Matrix The formula for calculating Mahalanobis distance is:

[0156]

[0157] in The larger the Markov distance, the further the current state deviates from the healthy domain, and the more severe the degradation of the shock resistance.

[0158] In this embodiment, to improve the robustness of the health domain boundary, a robust covariance estimation method is employed when establishing the health domain boundary. Specifically, the minimum covariance determinant (MCD) method is used to estimate the covariance matrix of the health data. The minimum covariance determinant method can effectively suppress the influence of a small number of outliers in the health data, thereby improving the accuracy of the health domain boundary.

[0159] For the 151st trajectory point, extract the previous point. and the next point (This point is collected during the subsequent start-up and shutdown process.)

[0160] Calculate vectors Vector magnitude Calculate vectors Vector magnitude .

[0161] Calculate the cosine of the angle between the vectors:

[0162]

[0163] Obtain the included angle radian.

[0164] Calculate local curvature .

[0165] The first 100 start-stop cycles of the spindle represent the healthy phase data, and the average of high-frequency abnormal indicators is statistically analyzed. Standard deviation Preset threshold During the 8-hour work shift on November 15, 2024, a total of 52 start-ups and shutdowns were executed, of which 7 were due to high-frequency abnormal indicators exceeding the threshold. The frequency of mutation events was calculated. The average frequency of mutation events in a healthy state is 0.35 times per hour. The current frequency is 2.5 times that of a healthy state, indicating an accelerated degradation phase.

[0166] Establish the health domain boundary based on the first 100 health status data. Calculate the mean vector of the three-dimensional joint feature space. The covariance matrix is:

[0167]

[0168] For the current feature vector Calculate the deviation vector .

[0169] Calculate the inverse of the covariance matrix to obtain the Mahalanobis distance:

[0170]

[0171] For the saliency level in three-dimensional space The critical value of the chi-square distribution is The radius of the health domain boundary is The current Mahalanobis distance of 8.67 exceeds the boundary of the healthy domain, indicating that the principal axis is in a state of significant degradation.

[0172] Table 6. Local characteristic analysis of the joint feature space trajectory:

[0173]

[0174] Step 800: Predict the growth trend of the topology sensitivity gradient using the damage accumulation criterion and crack propagation formula, calculate the remaining number of start-stop cycles before reaching the critical threshold, and generate the remaining lifetime prediction results.

[0175] The cumulative fatigue damage is calculated using the Palmgren-Miner linear damage accumulation criterion, based on the historical equivalent impact stress amplitude sequence calculated in step 600. The calculation formula is as follows:

[0176]

[0177] in, For stress amplitude category index, This represents the total number of stress amplitude categories. For the first The actual number of stress amplitude cycles, The fatigue life at this stress amplitude is determined by the SN curve.

[0178] By combining the Paris crack propagation formula, the growth trend of the norm of the topological sensitivity gradient tensor is predicted. The Paris formula is:

[0179]

[0180] in, Let be the norm of the topological sensitivity gradient tensor. For the number of start-stop cycles, The range of equivalent impact stress amplitude. and is a material constant.

[0181] Furthermore, the material constants in the aforementioned Paris formula... and The method involved fitting accelerated degradation test data of the same model of spindle to determine the optimal value. Specifically, the variation of the topological sensitivity gradient tensor norm with the number of start-stop cycles was measured under different equivalent impact stress amplitudes. Logarithmic linear regression was then performed on the experimental data to obtain the optimal value. and For the ceramic hybrid bearing spindle in this embodiment, the material constant... The typical range of values ​​is to ,index The typical value range is 2 to 4.

[0182] Based on the current cumulative fatigue damage and topological sensitivity gradient growth trend, the remaining number of start-stop cycles before reaching the critical sensitivity threshold or cumulative damage threshold is calculated, and the remaining lifetime prediction results are generated.

[0183] Furthermore, the specific steps for calculating the remaining number of start-stop cycles are as follows: First, numerically integrate the norm of the topological sensitivity gradient tensor using the Paris formula, and then convert the differential equation... Rewritten in integral form ,in Let the norm of the current topological sensitivity gradient tensor be denoted as . The critical sensitivity threshold, For the remaining number of start-stop cycles, assume the equivalent impact stress amplitude range. In future start-up and shutdown processes, the number of remaining start-up and shutdown cycles is kept at the average of the most recent 10 cycles. The remaining number of start-up and shutdown cycles is obtained by integrating the results and calculating the number of cycles based on topology sensitivity. Simultaneously, the remaining number of start-stop cycles based on cumulative damage is calculated using the Palmgren-Miner criterion. ,in The current cumulative fatigue damage, The average fatigue life under the current stress amplitude is used as the average fatigue life, and the smaller of the two values ​​is taken as the conservative prediction result for the remaining start-stop cycles. .

[0184] Furthermore, the aforementioned average fatigue life The calculation method is as follows: the equivalent impact stress amplitude of the most recent 10 start-stop processes is averaged to obtain the current average stress amplitude. Through the empirical formula of the SN curve Calculate the average fatigue life under this stress amplitude, where and The material parameters for the SN curve are: mean fatigue life, which represents the total number of start-stop cycles that the spindle can theoretically withstand under the current stress level.

[0185] Furthermore, the aforementioned critical sensitivity threshold Based on historical failure case statistics, specifically, topological sensitivity gradient data of at least 10 identical spindles were collected throughout their entire lifecycle from healthy to failed states. The gradient tensor norm before failure was extracted as a critical value. The mean of these critical values ​​minus one standard deviation was taken as a conservative critical sensitivity threshold, and the cumulative damage threshold was determined. According to the Palmgren-Miner criterion, the value is 1.0. When the cumulative damage reaches 1.0, it means that the spindle has theoretically reached its fatigue life limit.

[0186] It should be noted that the determination of the aforementioned SN curve refers to establishing the relationship between stress amplitude and fatigue life by fitting fatigue test data of spindles of the same model. The empirical formula for the SN curve is as follows: ,in and These are material parameters.

[0187] Furthermore, the material parameters of the aforementioned SN curves and The method for determining the fatigue life of spindles of the same model under different stress amplitudes is as follows: fatigue life tests are conducted on spindles of the same model under different stress amplitudes, and the number of failure cycles under each stress amplitude is recorded. and After taking the logarithm, perform linear regression to fit the parameters. and For ceramic hybrid bearing spindles, parameters The typical range of values ​​is to ,parameter The typical value range is 3 to 5.

[0188] In this embodiment of the application, to improve the accuracy of remaining lifetime prediction, the effect of stress ratio is considered when applying the Paris formula. The modified Paris formula is:

[0189]

[0190] in, Stress ratio, and These are the minimum and maximum values ​​of the stress cycles, respectively. This is the stress ratio influence index.

[0191] Furthermore, the aforementioned stress ratio influence index The method was determined by fitting fatigue test data, specifically: under different stress ratios Fatigue tests were conducted under specific conditions to measure the relationship between crack propagation rate and stress ratio. The experimental data were then fitted to obtain an exponential coefficient. For ceramic hybrid bearing spindles, the stress ratio influence index The typical value range is 0.5 to 1.5. The stress ratio influence index reflects the degree of influence of the average stress level of the stress cycle on the crack propagation rate.

[0192] Rainflow counting analysis was performed on the time series of thermal-electric phase space trajectory curvature extracted during the start-up and shutdown process. The curvature series contains 28 sampling points, with a maximum curvature of 3.67 and a minimum curvature of 0.15. The rainflow counting algorithm was applied to extract the extreme value sequence, identifying seven complete stress cycles. The amplitude of the first type of cycle was... It appears once; second type of cyclic amplitude It appears twice; the third type of cyclic amplitude It appeared 4 times. Total number of loops. The reference trajectory length under the healthy state of the spindle. Unit length, current trajectory length Influence coefficient of unit length and trajectory length Reference duration seconds, current transient duration Seconds, duration affects index The equivalent impact stress amplitude is calculated using the formula:

[0193]

[0194] This spindle has accumulated 150 start-stop data points, recorded from the initial stage of operation. The Frobenius norm of the topology sensitivity gradient tensor is used. , condition number and equivalent impact stress amplitude A joint feature space is constructed by combining three features. Mean normalization based on the range is then applied to each of the three features: (Based on 150 historical start-stop cycles). The range is from 2.15 to 5.82.

[0195] The current value is 3.94, normalized to ; The range is from 2.38 to 8.94, and the current value of 5.68 is normalized to... ; The range is from 1.45 to 3.67, and the current value of 2.24 is normalized to... Constructing a 3D joint feature vector This vector becomes the 151st point of the joint feature space trajectory.

[0196] Table 7. Rainflow Count Cycle Classification Statistics:

[0197]

[0198] Table 8. Joint feature space construction data:

[0199]

[0200] Figure 5 The study shows the evolution trend of cumulative fatigue damage of the spindle from the start of operation to the 151st start-stop process.

[0201] Figure 6 The evolution trends of the topology sensitivity gradient tensor norm G (left Y-axis, broken line) and the predicted number of remaining starts and stops N_remain (right Y-axis, bar chart) are shown during the process of the main axis starting and stopping from the 100th to the 155th time.

[0202] Step 900: Determine whether the remaining number of start-stop cycles meets the production plan requirements, and generate a warning signal for the degradation of the spindle's shock resistance, a ceramic bearing replacement suggestion, and a start-stop frequency optimization strategy;

[0203] Obtain the production plan from the CNC system and extract the expected number of start-stop cycles within the planned period. Determine if the remaining number of start-stop cycles calculated in step 800 is greater than or equal to the expected number of start-stop cycles. If the remaining number of start-stop cycles is less than the expected number of start-stop cycles, generate a spindle shock resistance degradation warning signal and output it to the CNC system's human-machine interface and maintenance management system. Simultaneously, based on the analysis results of the current topology sensitivity gradient tensor, if the high-frequency band topology sensitivity gradient increases significantly, generate a ceramic bearing replacement suggestion. If the remaining number of start-stop cycles is close to the expected requirement, generate a start-stop frequency optimization strategy, suggesting reducing the number of start-stop cycles or extending the start-stop interval through process optimization, and output it to the production scheduling system.

[0204] Furthermore, the aforementioned criterion for a significant increase in the high-frequency band topology sensitivity gradient refers to the current high-frequency band topology sensitivity gradient. The norm exceeds the mean of the gradient norm in a healthy state. Add 2 standard deviations ,Right now The judgment criterion for a significant increase in the high-frequency band topological sensitivity gradient is based on the 2σ criterion, indicating that the current high-frequency band gradient is in an abnormal region of healthy distribution, corresponding to increased microscopic damage on the surface of the ceramic bearing.

[0205] It should be noted that the above-mentioned start / stop frequency optimization strategy is generated by calculating the recommended reduction percentage of start / stop times based on the ratio of remaining start / stop times to expected demand. If the remaining start / stop times are 50%-80% of expected demand, a reduction of 20% is recommended; if the remaining start / stop times are 30%-50% of expected demand, a reduction of 40% is recommended; if the remaining start / stop times are less than 30% of expected demand, maintenance should be scheduled immediately.

[0206] Furthermore, the percentage thresholds in the aforementioned start-stop frequency optimization strategy are determined through production practice experience and risk assessment. When the remaining start-stop counts are 50%-80% of the expected demand, it indicates that the spindle is in a mild degradation stage. Reducing the start-stop counts by 20% can slow down the degradation process while ensuring production efficiency. When the remaining start-stop counts are 30%-50% of the expected demand, it indicates that the spindle is in a moderate degradation stage. It is necessary to reduce the start-stop counts by 40% to avoid accelerating damage. When the remaining start-stop counts are less than 30% of the expected demand, it indicates that the spindle is in a severe degradation stage. Continuing to operate poses a risk of sudden failure, and maintenance should be arranged immediately.

[0207] In this embodiment, to achieve closed-loop control for spindle fault prediction, a spindle protection control strategy is automatically triggered after the early warning signal is generated. Specifically, when severe degradation of shock resistance is detected (remaining start-stop cycles are less than 20% of the expected requirements), the maximum spindle speed and maximum cutting load are automatically limited, restricting the maximum spindle speed to 70% of the rated speed and the maximum cutting power to 60% of the rated power, to prevent the spindle from being subjected to excessive load and sudden failure in a degraded state.

[0208] Furthermore, the speed and power limiting ratios in the aforementioned spindle protection control strategy are determined through safety margin analysis. Limiting the maximum speed to 70% of the rated speed reduces the contact stress of the ceramic bearing to a safe range, while limiting the maximum cutting power to 60% of the rated power keeps the radial and axial loads on the spindle below the safety threshold under degradation conditions. These limiting ratios comprehensively consider the remaining load-bearing capacity of the spindle under degradation conditions and the basic requirements of production and processing, ensuring that the spindle does not experience sudden failure while retaining a certain processing capacity to complete emergency production tasks.

[0209] The production plan for the spindle for the next 30 days is obtained from the CNC system. During the planned period, four batches of aluminum alloy parts need to be processed, with an average processing cycle of 6 days per batch. The daily working hours are 8 hours, and each shift has an average of 50 start-stop operations. The expected total number of start-stop operations is [missing information]. The predicted remaining number of start-stop cycles (13,770) is compared with the expected demand (1,200), and the ratio is calculated as follows: The value is significantly greater than 1, indicating that the spindle can operate normally within the planned period. However, considering that the current Mahalanobis distance of 8.67 exceeds the healthy domain boundary radius of 3.37, and the mutation event frequency of 0.875 times / hour is 2.5 times that of the healthy state, the system determines that the spindle is in a moderately accelerated degradation phase and generates preventative maintenance recommendations.

[0210] The high-frequency band topological sensitivity gradient was examined, and the maximum component of the high-frequency band differential gradient vector calculated in step 400 was found to be 1.68. The mean norm of the high-frequency band gradient under healthy conditions was statistically analyzed. Standard deviation The threshold for judgment is The current maximum gradient component in the high-frequency band is 1.68, which exceeds the threshold of 1.56 and is considered to be significantly increased. The system automatically generates a ceramic bearing replacement recommendation, suggesting that the bearing replacement be performed within the maintenance window after the completion of the current batch of processing.

[0211] Based on the remaining start / stop cycles and the current degradation state, the system generates a start / stop frequency optimization strategy. Although there are sufficient remaining start / stop cycles, to slow down the degradation process, the system recommends optimizing process parameters to reduce unnecessary start / stop operations. By analyzing historical machining data, it was identified that approximately 15% of the start / stop operations per shift (about 7-8 times) were due to short-term shutdowns followed by restarts caused by unreasonable process switching. It is recommended to reduce the number of start / stop cycles per shift from 50 to 42-43 by optimizing the workpiece clamping sequence and tool pre-setting strategy, a reduction of approximately 15%, which is expected to reduce the spindle degradation rate by 12%.

[0212] The system outputs early warning information and optimization suggestions to the CNC system's human-machine interface, displaying "Moderate degradation warning for spindle impact resistance; it is recommended to replace the ceramic bearing within the 6-day maintenance window. The current remaining lifespan is 13,770 start-stop cycles, which can meet recent production needs." Simultaneously, detailed degradation index data is sent to the maintenance management system, including key parameters such as Mahalanobis distance (8.67), abrupt event frequency (0.875 times / hour), and cumulative damage degree (0.00574), for maintenance personnel to develop maintenance plans. Start-stop frequency optimization strategies are sent to the production scheduling system, suggesting adjustments to process parameters to reduce start-stop frequency.

[0213] Table 9. Results of Early Warning Decision-Making and Strategy Generation:

[0214]

[0215] This implementation decomposes the current envelope into three frequency band subsequences—low-frequency, mid-frequency, and high-frequency—using a multi-resolution filter bank. It independently calculates the continuous homology topology features for each frequency band, overcoming the factor that low-frequency components mask weak high-frequency signals in full-band topology analysis. This allows it to capture subtle topological changes caused by high-frequency degradation symptoms such as microcracks in ceramic bearings.

[0216] This implementation calculates the differential gradient of the topological feature vectors of each frequency band before and after the transient, generates a topological sensitivity gradient tensor, and quantifies the dynamic response characteristics of the spindle system topology during the transient process. This overcomes the factor that steady-state topology monitoring ignores the dynamic response during the transient process, thereby enabling the identification of the degradation of the spindle's topological recovery capability under start-stop transient impacts.

[0217] This implementation method constructs a two-dimensional phase space trajectory of thermo-electricity and uses a rainflow counting algorithm to quantify the equivalent impact stress amplitude of each start-up and shutdown, thereby achieving a quantitative assessment of the transient thermo-electric impact intensity and overcoming the limitations of independently analyzing temperature or current to comprehensively reflect the impact intensity.

[0218] This implementation constructs a joint feature space of the topological sensitivity gradient tensor and the equivalent impact stress amplitude, and calculates the Mahalanobis distance from the joint space trajectory to the boundary of the healthy domain. This enables the coupled analysis of the dynamic response characteristics of the topological structure and the transient impact intensity, overcoming the limitation that independent analysis of topological features or impact intensity cannot reflect the coupling effect between the two. Thus, it can accurately assess the sensitivity changes of the spindle system to impact.

[0219] This implementation utilizes the Palmgren-Miner linear damage accumulation criterion to quantify the cumulative fatigue damage during historical start-stop processes, and combines it with the Paris crack propagation formula to predict the growth trend of the topological sensitivity gradient. This enables quantitative modeling of the spindle's shock resistance degradation process and prediction of its remaining life, thereby providing accurate early warning before the spindle loses its shock resistance under the combined effects of multi-band topological degradation and cumulative fatigue damage, and preventing sudden failure of the spindle under frequent start-stop conditions.

[0220] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A numerical control machining center spindle fault prediction control method, characterized in that, The method comprises the following steps: Obtaining current data and temperature data of the main shaft motor during the start-stop process, identifying the data segment type based on the start-stop instruction timestamp, and generating a labeled data set; Performing Hilbert transform on the three-phase current data of each stage and extracting the current envelope line, decomposing the envelope line into a plurality of frequency band subsequences using a multi-resolution filter bank, and generating a frequency band subsequence set; Constructing a topological space for each frequency band subsequence and calculating persistent homology features to generate a multi-frequency band topological feature matrix; Calculating a high-frequency anomaly index and a topological sensitivity gradient tensor based on the multi-frequency band topological feature matrix to generate a topological dynamic response feature; Constructing a thermal-electric two-dimensional phase space trajectory based on the temperature rise rate and the current amplitude coefficient of variation, calculating the geometric features of the trajectory, and generating a transient impact feature; Quantifying the transient impact strength using the rain flow counting algorithm, constructing a joint feature space of the topological sensitivity gradient tensor and the equivalent impact stress amplitude, and generating a joint feature space trajectory; Calculating the curvature of the joint feature space trajectory and the mutation event frequency of the high-frequency anomaly index, calculating the Mahalanobis distance based on the health domain boundary, and generating an impact resistance degradation index; Using the damage accumulation criterion and the crack propagation formula to predict the growth trend of the topological sensitivity gradient, calculating the remaining number of start-stops to reach the critical threshold, and generating a remaining life prediction result.

2. The numerical control machining center spindle fault pre-judgment control method according to claim 1, characterized in that, The step of constructing a topological space for each frequency band subsequence and calculating persistent homology features comprises: For high-frequency sub-sequences, a Takens delay embedding theorem is used to construct a delay embedding phase space, and a one-dimensional time series is reconstructed into a point set in an m-dimensional phase space; The optimal delay time is determined by the mutual information method, and the delay value corresponding to the first local minimum value of the mutual information curve is selected; The optimal embedding dimension is determined by the pseudo-neighbor method, and when the pseudo-neighbor ratio falls below 5%, the corresponding embedding dimension is the optimal value; The point set of the phase space is topologically connected at different scales, and a Vietoris-Rips complex construction method is used to form a simple complex sequence; The generation and disappearance of topological features at each scale are calculated to generate a persistent homology barcode; The statistical features of the barcode length are extracted as topological invariants, including mean, variance, maximum, and entropy, to form a topological feature vector.

3. The numerical control machining center spindle fault pre-judgment control method according to claim 1, characterized in that, In the step of calculating a high-frequency anomaly index based on a multi-frequency band topological feature matrix, the calculation formula of the high-frequency anomaly index is the energy ratio of the high-frequency topological feature vector to the energy of the medium-frequency and low-frequency topological feature vectors, that is, the square of the Euclidean norm of the high-frequency topological feature vector divided by the sum of the squares of the Euclidean norms of the medium-frequency and low-frequency topological feature vectors.

4. The numerical control machining center spindle fault pre-judgment control method according to claim 1, characterized in that, The topological sensitivity gradient tensor is a difference gradient matrix between the topological feature vectors of each frequency band before and after the transient state, which contains the difference values of the topological feature vectors of the low-frequency band, the medium-frequency band, and the high-frequency band. The norm of the gradient tensor reflects the overall amplitude of the change in the topological structure, and the condition number reflects the uneven degree of topological changes between different frequency bands.

5. The numerical control machining center spindle failure pre-judgment control method according to claim 1, characterized in that, The step of constructing a thermal-electric two-dimensional phase space trajectory based on the temperature rise rate and the current amplitude coefficient of variation comprises: The first derivative of the temperature time series is calculated using the central difference method as the temperature rise rate; The ratio of the standard deviation to the mean of the current amplitude time series is calculated as the coefficient of variation; The temperature rise rate and the current variation coefficient are respectively subjected to Z-score standardization processing; A two-dimensional phase space trajectory is constructed based on the normalized thermal state variable and the electrical state variable; The maximum curvature and the trajectory length of the phase space trajectory are calculated as geometric features.

6. The numerical control machining center spindle fault pre-judgment control method according to claim 5, characterized in that, The calculation method of the maximum curvature of the trajectory is as follows: the first and second derivatives of the thermal state variable and the electrical state variable of the phase space trajectory are calculated, the curvature values at each time are obtained according to the curvature calculation formula, and the maximum value is taken as the maximum curvature of the trajectory; when the trajectory speed is less than a preset tolerance, the curvature at this time is set to zero.

7. The numerical control machining center spindle failure pre-judgment control method according to claim 1, characterized in that, The step of quantifying the transient impact strength by using the rainflow counting algorithm comprises: The curvature time sequence of the thermal-electric coupling phase space trajectory is taken as an equivalent stress time history; The peak and valley values of the curvature time sequence are extracted to form an extreme sequence; Complete cycles are identified and the curvature amplitude of each cycle is calculated; The number of cycles in different amplitude ranges is counted; The equivalent impact stress amplitude of a single start-stop is calculated by weighting the cycle amplitude in combination with the trajectory length and the transient duration.

8. The numerical control machining center spindle failure pre-judgment control method according to claim 1, characterized in that, The step of calculating the Mahalanobis distance based on the health domain boundary comprises: At least 100 start-stop processes of the main shaft health state are collected to obtain joint feature data; The mean vector and the covariance matrix of the health data are calculated; Robust covariance estimation is performed by using the minimum covariance determinant method; The health domain boundary is defined as a hyperellipsoid with the Mahalanobis distance equal to the critical value of the chi-square distribution; The Mahalanobis distance of the current joint feature vector relative to the mean and the covariance matrix of the health domain is calculated.

9. The numerical control machining center spindle failure pre-judgment control method according to claim 1, characterized in that, The step of predicting the growth trend of the topology sensitivity gradient, calculating the remaining start-stop number reaching the critical threshold value by using the damage accumulation criterion and the crack propagation formula comprises: The cumulative fatigue damage degree of the historical start-stop process is calculated based on the Palmgren-Miner linear damage accumulation criterion; The growth trend of the norm of the topology sensitivity gradient tensor with the number of start-stops is predicted by using the Paris crack propagation formula; The remaining start-stop number for the topology sensitivity gradient to reach the critical threshold value is calculated by numerical integration; The average fatigue life under the current stress level and the remaining start-stop number based on the cumulative damage are calculated based on the S-N curve; The smaller value of the topology sensitivity prediction and the cumulative damage prediction is taken as the final prediction result of the remaining start-stop number.

10. A numerical control machining center spindle fault prediction control system for performing the numerical control machining center spindle fault prediction control method according to any one of claims 1-9, characterized in that, It comprises: A data acquisition module for acquiring current data and temperature data of the main shaft motor start-stop process; A multi-band decomposition module for multi-resolution filter decomposition of the current envelope; A topology feature extraction module for calculating the persistent homology topology features of each frequency band; A dynamic response analysis module for calculating the topology sensitivity gradient tensor and the high-frequency anomaly index; A thermal-electric coupling analysis module for constructing a thermal-electric two-dimensional phase space trajectory and quantifying the impact strength; A joint feature analysis module for constructing a joint feature space and calculating the Mahalanobis distance; A life prediction module for predicting the remaining start-stop number based on damage accumulation and crack propagation; A decision output module for generating a warning signal and an optimization strategy.

Images