A method for measuring dynamic and static contact stiffness of bearing rollers based on digital twin model

Abstract

The present invention discloses a method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model, comprising: separately constructing a rotor system and a bearing system; for the rotor system, designing an algorithm based on VGWO‑DBSCAN‑AMD to perform modal parameter identification and designing a response surface algorithm to perform simulation model correction, thereby obtaining a high-precision rotor system simulation model; for the bearing system, conducting a friction test to obtain the bearing roller friction coefficient and building an initial simulation model of the bearing system; obtaining a digital twin model by assembling the bearing system and rotor system simulation models; operating a system test bench and the digital twin model under the same operating conditions and simultaneously extracting time domain vibration signals, performing correlation analysis on the vibration characteristics using a correlation evaluation criterion, adjusting the roller contact stiffness to obtain the true value of the contact stiffness of the test bench under different operating conditions, and fitting the contact stiffness under multiple operating conditions to obtain the bearing roller contact stiffness value within the full operating range of the test bench.

Description

Technical Field

[0001] The present invention relates to the field of bearing roller dynamic and static contact stiffness measurement, and in particular to a bearing roller dynamic and static contact stiffness measurement method based on a digital twin model. Background Art

[0002] The dynamic and static contact stiffness of bearing rollers is a crucial indicator of bearing operating conditions. Obtaining accurate contact stiffness is crucial for studying bearing motion characteristics, condition monitoring, and health maintenance. Traditional measurement methods primarily focus on static contact stiffness, with limited research on dynamic contact. In reality, the dynamic contact stiffness of bearing rollers varies with changes in equipment load and speed. Because rotating equipment cannot be directly equipped with measuring devices for stiffness testing, a mature method for measuring dynamic bearing stiffness currently exists.

[0003] Digital twin technology can be used to measure the dynamic and static contact stiffness of bearings using vibration signals. While ensuring the accuracy of all parameters except the dynamic and static contact stiffness parameters of the bearing rollers, the dynamic and static contact stiffness of the bearing rollers can be adjusted to achieve similar time domain vibration signal characteristics of the simulation and test systems under the same working conditions, thereby obtaining accurate dynamic and static contact stiffness of the bearings. Summary of the Invention

[0004] In view of this, the purpose of the present invention is to provide a method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model, aiming to measure the dynamic and static contact stiffness of bearings under complex and variable working conditions within the entire operating range of the equipment based on digital twin technology, and to improve the generalization performance, accuracy, robustness and reliability of the measurement method.

[0005] In order to achieve the above object, the present invention adopts the following technical solutions:

[0006] A method for measuring dynamic and static contact stiffness of bearing rollers based on a digital twin model, the method comprising:

[0007] Step S1: constructing a bearing system geometry simulation model and a rotor system geometry simulation model for a rotor-bearing system test bench, and then performing a hierarchical modal test on each component in the rotor system geometry simulation model;

[0008] Step S2: Based on the results of the hierarchical modal test in step S1, the modal parameters are identified using the VGWO-DBSCAN-AMD algorithm, and the modal parameter values ​​are input into the ADAMS finite element simulation model to construct a rotor system simulation model; wherein the method includes: firstly determining the optimal modal number based on the improved DBSCAN clustering algorithm, then performing noise reduction processing on the vibration signal based on the optimal modal number and the VMD algorithm, then determining the boundary splitting frequency using the autoregressive power spectrum, then determining the optimal splitting frequency using VGWO and performing modal decomposition based on AMD to obtain modal components of each order, and finally obtaining the modal vibration shape;

[0009] Step S3: For the rotor system simulation model obtained in step S2, a response surface model correction algorithm based on VGWO-DBSCAN-SVM is used to correct it, and a high-precision rotor system simulation model is constructed according to the corrected model parameters;

[0010] Step S4: performing friction experiments and initializing the dynamic and static contact stiffness of the bearing rollers on the bearing system geometry simulation model to build an initial simulation model of the bearing system, wherein the dynamic and static friction coefficients of the bearings are obtained through the friction experiments;

[0011] Step S5: constructing a digital twin model of the rotor-bearing system based on the high-precision rotor system simulation model obtained in step S3 and the initial bearing system simulation model obtained in step S4;

[0012] Step S6: Run the rotor-bearing system test bench and the digital twin model established in step S5 under the same complex variable operating conditions, extract the output vibration signal characteristics, and perform similarity evaluation using a correlation evaluation criterion;

[0013] Step S7: For the similarity evaluation result obtained in step S6, if the similarity of the vibration signal characteristics is lower than a preset value, the dynamic and static contact stiffness of the bearing roller in the simulation model is adjusted until the similarity of the output vibration signal characteristics is higher than the preset value, and the current bearing roller contact stiffness is output. The bearing roller contact stiffness value within the full operating range of the test bench is obtained by fitting the bearing roller contact stiffness values ​​under multiple working conditions.

[0014] Furthermore, the step S1 includes:

[0015] Step S101, constructing a bearing system geometric simulation model, which includes: establishing a bearing system geometric simulation model based on ADAMS for the bearing cage, rollers, and inner and outer rings of the rotor-bearing system test bench;

[0016] Step S102, constructing a rotor system geometric simulation model, which includes: constructing an ADAMS-based rotor system geometric simulation model by removing components other than the components used to construct the ADAMS-based bearing system geometric simulation model for the rotor-bearing system test bench;

[0017] Step S103: performing modal tests on each component in the rotor system geometric simulation model and the rotor system itself.

[0018] Furthermore, in step S2, determining the optimal modal number based on the improved DBSCAN clustering algorithm specifically includes:

[0019] Step S2011: Calculate the vibration response data set obtained based on the layered modal test in step S1 The local density of the data points in is as follows:

[0020]

[0021] In formula (1), d c represents the cutoff distance; η is the number of data sets; I S ={1,2,...,η} is the indicator set corresponding to the data set; d ij For data point γ i and γ j The distance between i For S and γ i The distance between them is less than d c The number of data points;

[0022] Step S2012: Design express A descending subscript order of q1 ≥ρ q2 ≥…≥ρ qη ; Calculate the distance as follows:

[0023]

[0024] In formula (2), Represents data points With data points The distance between When the local density is the largest, δ qi Indicates that the data set is The data point with the largest distance The distance between them; in all local densities greater than Among the data points, The data point with the smallest distance the distance between them;

[0025] Step S2013, determining the cluster center, which includes: selecting points with high local density and relatively high distance as the cluster center, assigning the remaining points to the corresponding nearest neighbor cluster with higher density, and determining the optimal mode number K=3.

[0026] Furthermore, in step S2, the noise reduction process is performed on the vibration signal based on the optimal modal number and the VMD algorithm, which specifically includes:

[0027] Decomposing and denoising the vibration signal extracted in step S1 according to the optimal modal number K obtained in step S2013;

[0028] The Pearson correlation between the components after VMD decomposition and the original signal is analyzed as the basis for denoising and reconstruction.

[0029] Furthermore, in step S2, the use of the autoregressive power spectrum to determine the boundary segmentation frequency specifically includes:

[0030] The AR model spectrum is used instead of the Fourier spectrum to determine the power spectrum peak and the boundary division frequency of AMD. The AR model represents the signal as follows:

[0031]

[0032] In formula (3), p is the order of the AR model; u(n) has a mean of 0 and a variance of White noise sequence; a k (k=1,2,...,p) is the corresponding p-order model parameter;

[0033] The standard equation of the autoregressive model is as follows:

[0034]

[0035] In formula (4), r x (m) is the autoregressive variable; r k (k) is the autoregressive variable sequence;

[0036] Then the autoregressive power spectrum P of the signal x(n) is AR (e jω ) can be calculated by the following formula:

[0037]

[0038] The AR model spectrum is estimated using the Burg method.

[0039] Furthermore, in step S2, determining the optimal splitting frequency using VGWO and performing modal decomposition based on AMD specifically includes:

[0040] Step S2041: Determine the first D-order damped natural frequencies of the structure ω=[ω1, ω2, ..., ω D ] T ;

[0041] Step S2042: Use VGWO to optimize the boundary segmentation frequency in AMD to obtain the optimal vibration attenuation curve, and create a population D consisting of n particles based on VGWO = (D1, D2, ..., D n ), where D i =[d i1 ,d i2 ,...,d iD ] T , then the boundary division frequency of AMD can be expressed as:

[0042]

[0043] In formula (6), d ij is the boundary frequency cutoff bandwidth of the j-th mode of the structure in the i-th particle; n is the number of particles;

[0044] Step S2043: In order to optimize the parameters using the performance index, the design fitness value can be obtained by decomposing and separating the signal x from AMD. j (t) and the remaining signal x k The correlation coefficient of (t) is calculated and the fitness function satisfies:

[0045]

[0046] In formula (7), ρ j is x j (t) and x k Pearson correlation coefficient of (t); and They are samples x j (t) and x k (t) mean; T is the number of samples; if x j (t) is completely separated, then x j (t) and x k The correlation coefficient of (t) is the smallest;

[0047] Step S2044: Introduce the speed component of the particle swarm optimization algorithm into the gray wolf optimization algorithm to form a variable speed gray wolf optimization algorithm. The speed and position components of the fused particle swarm algorithm are as follows:

[0048]

[0049] p i (m+1)=pi (m)+v i (m+1) (9)

[0050] In formulas (8) and (9), v i is the speed of the i-th gray wolf; p i is the current position of the i-th gray wolf; c1, c2, c3 are learning factors and satisfy c1, c2, c3∈[0,1]; ζ is the inertia factor; m is the number of iterations; X GWO1 、X GWO2 、X GWO3 Gray Wolf GWO Relative to Gray Wolf α GWO , β GWO , δ GWO The forward vector of

[0051] The obtaining of the modal vibration shape specifically includes:

[0052] After decomposing the modal components of each order by AMD in step S204, the modal vibration shape is obtained by considering the structural modal response at all l sensor positions, where the k-th modal vibration shape vector of the structure is as follows:

[0053]

[0054] In formula (10), u k,j (t m ) is the kth modal response at the jth sensor position; t m is the time when the local minimum or local maximum of the modal response occurs; l is the number of installed vibration sensors; Nor(·) represents the normalized calculation; take k = 1, 2, ..., K respectively, and repeat the calculation of formula (10) to obtain the modal vibration matrix

[0055] Furthermore, the step S3 specifically includes:

[0056] Step S301: determining the optimal modal number K based on the improved DBSCAN clustering algorithm;

[0057] Step S302: Use support vector machine to fit the response surface to establish an SVM response surface model between the correction variable Γ and the modal frequency. The number of response surface fittings is set to K according to the optimal mode number obtained in step S301. Given l sample sets of correction variable Γ and modal frequency:

[0058] {(Γ1,Θ 1,k ),(Γ2,Θ 2,k ),...,(Γ l ,Θ l,k )} (11)

[0059] In formula (11), l is the number of samples; Γ i is the variable to be corrected; Θ i,k is the sample output; k = 1, 2, ..., K is the k-th modal frequency; with the fitting accuracy ε as the error function, the optimization problem is as follows:

[0060]

[0061] In formula (12), α i and is the Lagrange coefficient; κ(Γ i ,Γ j ) is any kernel function that satisfies the Mercer condition;

[0062] Solving equation (12) yields the nonlinear regression fitting function of the kth-order modal frequency and the correction variable Γ:

[0063]

[0064] Step S303: using the VGWO algorithm to perform model correction, which specifically includes:

[0065] Based on the support vector machine response surface established in step S302, the measured modal frequency and the theoretical modal frequency Θ i,k The objective function is established as follows:

[0066]

[0067] In formula (14), Γ L and Γ U are the upper and lower limits of the correction variable, respectively.

[0068] Furthermore, the step S5 specifically includes:

[0069] By using ADAMS simulation software, the outer ring and inner ring of the bearing in the initial simulation model of the bearing system obtained in step S4 are respectively connected to the bearing seat in the high-precision rotor system simulation model obtained based on step S3 and to the rotating shaft through external connection points; at the same time, the bearing rollers are set to contact with the inner and outer rings, a revolute pair is added to the rotating shaft, and a thin layer unit is added to the platform base and fixedly connected to the ground.

[0070] Furthermore, the step S6 specifically includes:

[0071] Step S601: performing vibration signal feature extraction, which includes:

[0072] Step S6011: For the rotor-bearing system test bench and the digital twin model, the load and speed of the rotor-bearing system are simultaneously varied. The load and speed variation range should cover the operating range of the rotor-bearing system test bench. Time domain vibration signals under complex variable operating conditions are extracted at the same location on the test bench and the digital twin model using the same sampling frequency.

[0073] Step S6012: extracting the following features from the acquired variable working condition vibration data: time domain features, frequency domain features, and time-frequency domain features;

[0074] Step S6013: construct a feature set using the above features;

[0075] Step S602: Execute feature similarity evaluation, which specifically includes:

[0076] The time domain features, frequency domain features, and time-frequency domain features of the test bench and the digital twin model obtained in step S601 are similarly evaluated using the correlation evaluation criteria. The correlation evaluation criteria are constructed as follows:

[0077]

[0078] In formula (15), F h and l h Represent the characteristic value of the h-th sample and the corresponding time sequence respectively; and are the sample eigenvalue sequence and the time series mean, respectively; H is the number of samples; the value of the correlation evaluation index ranges from 0 to 1. The better the characteristic correlation between the test bench and the digital twin model vibration signal, the closer the value is to 1, otherwise it is closer to 0.

[0079] Furthermore, the step S7 specifically includes:

[0080] Step S701: Obtain the dynamic and static contact stiffness of the bearing roller under a single working condition, which specifically includes:

[0081] When the similarities of all features obtained in step S602 exceed 0.7, it indicates that the vibration signal characteristics of the test bench and the digital twin model are basically consistent, and the current dynamic and static contact stiffness of the bearing roller is considered to be the bearing roller contact stiffness under the corresponding working condition;

[0082] Otherwise, the dynamic and static contact stiffness of the bearing rollers is adjusted until all characteristic similarities of the output vibration signals exceed 0.7;

[0083] Step S702: acquiring the dynamic and static contact stiffness of the bearing rollers in the full operating range, which specifically includes:

[0084] The dynamic and static contact stiffness of the bearing rollers obtained under multiple working conditions with varying loads and speeds are fitted to obtain the bearing contact stiffness variation curve within the full operating range of the rotor-bearing test bench.

[0085] The beneficial effects of the present invention are:

[0086] 1. The present invention uses the VGWO-DBSCAN-AMD algorithm to identify modal parameters, effectively solving the problem of modal order determination and achieving high recognition accuracy for dense modes in strong interference environments;

[0087] 2. The response surface model correction based on VGWO-DBSCAN-SVM in the present invention can automatically identify the modal order, avoiding manual mode picking and significantly improving the accuracy of model correction;

[0088] 3. The present invention analyzes the vibration signal responses of the simulation and test systems based on feature similarity, avoiding the stringent requirement of complete consistency of the time domain vibration signals and having greater robustness and higher reliability against unknown disturbances and noise;

[0089] 4. The present invention measures the dynamic and static contact stiffness of bearing rollers using only the time-domain vibration response signal, which reduces the number and space required for the installation of measurement devices, and has a wider range of applications and higher generalization performance;

[0090] 5. The present invention measures the dynamic and static contact stiffness of the bearing rollers under various complex and variable working conditions of the equipment, and obtains a high-precision stiffness time-varying function within the full operating range of the equipment through curve fitting. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] Figure 1 This is a flow chart of a method for measuring the dynamic and static contact stiffness of a bearing roller based on a digital twin model provided in Example 1;

[0092] Figure 2 This is a schematic diagram of performing a layered modal test on a rotor system provided in Example 1, wherein: Figure 2 a is a field photo of the rotor system. Figure 2 b is the modal test structure diagram;

[0093] Figure 3 Schematic diagram of response data of the layered modal test provided in Example 1, wherein: Figure 3 a is a schematic diagram showing the magnitude of the knocking force at different sampling points. Figure 3 b represents a schematic diagram of the frequency response function;

[0094] Figure 4 Flowchart of the modal parameter algorithm based on VGWO-DBSCAN-AMD provided in Example 1;

[0095] Figure 5 This is a schematic diagram of the improved DBSCAN clustering results provided in Example 1;

[0096] Figure 6 Schematic diagram of the VMD noise reduction result provided in Example 1, Figure 6 a represents the original signal in the time domain, Figure 6 b represents the signal after VMD denoising and reconstruction;

[0097] Figure 7 Schematic diagram of the correlation analysis results between each component and the original signal after VMD decomposition provided in Example 1;

[0098] Figure 8 is a schematic diagram of the SVM-based fitting response surface provided in Example 1, wherein, Figure 8 a- Figure 8 The i's represent the response surfaces of the 1st to 9th order mode fittings respectively;

[0099] Figure 9 Schematic diagram of the model correction results provided in Example 1, wherein: Figure 9 a- Figure 9 c represents the model correction results of the 1st to 3rd order modes of the rotor system;

[0100] Figure 10 The measurement curve of the friction coefficient of the bearing roller in the embodiment of the present invention;

[0101] Figure 11 This is the time domain feature of the vibration signal provided in Example 1, wherein the meaning of each sub-graph is shown in the subscript of the figure;

[0102] Figure 12 This is the frequency domain feature of the vibration signal provided in Example 1, wherein the meaning of each sub-graph is shown in the subscript of the figure;

[0103] Figure 13 The time-frequency domain characteristics of the vibration signal provided in Example 1;

[0104] Figure 14 This is a schematic diagram of the dynamic and static contact stiffness of the bearing rollers in the full operating range of the equipment provided in Example 1. DETAILED DESCRIPTION

[0105] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0106] Example 1

[0107] See also Figures 1-14 This embodiment provides a method for measuring the dynamic and static contact stiffness of a bearing roller based on a digital twin model. The process of the method is as follows: Figure 1 As shown, the method specifically includes the following steps:

[0108] Step S1: constructing a bearing system geometry simulation model and a rotor system geometry simulation model for a rotor-bearing system test bench, and then performing a hierarchical modal test on each component in the rotor system geometry simulation model;

[0109] Specifically, in this embodiment, step S1 includes:

[0110] Step S101, constructing a bearing system geometric simulation model, which includes: establishing a bearing system geometric simulation model based on ADAMS for the bearing cage, rollers, and inner and outer rings of the rotor-bearing system test bench, wherein the basic structural parameters of the bearing system geometric simulation model are consistent with the basic structural parameters of the actual bearing in the test bench;

[0111] Step S102, constructing a rotor system geometric simulation model, which includes: constructing an ADAMS-based rotor system geometric simulation model for the rotor-bearing system test bench except for the components used to construct the ADAMS-based bearing system simulation model; more specifically, in this embodiment, the rotor system specifically includes: a platform base, a bracket, a bearing seat, a rotating shaft, an eccentric disk and other components.

[0112] Step S103: Perform modal tests on each component in the rotor system geometric simulation model and the rotor system itself. More specifically, the modal test model is as follows: Figure 2 The input and output data and frequency response function in the modal test are shown in Figure 3 As shown, in which Figure 3 In the layered modal test, a total of 12 knocking points were set for data collection.

[0113] Step S2: Based on the results of the hierarchical modal test in step S1, the modal parameters are identified using the VGWO-DBSCAN-AMD algorithm, and the modal parameter values ​​are input into the ADAMS finite element simulation model to construct a rotor system simulation model. The algorithm flow is as follows: Figure 4 As shown, the modal parameter values ​​are input into the ADAMS finite element simulation model;

[0114] Specifically, in this embodiment, step S2 includes:

[0115] Step S201: Determine the optimal modal number K based on the improved DBSCAN clustering algorithm, which includes:

[0116] Step S2011: Calculate the vibration response data set obtained based on the layered modal test in step S1 The local density of the data points in is as follows:

[0117]

[0118] In formula (1), d c represents the cutoff distance; η is the number of data sets; I S ={1,2,...,η} is the indicator set corresponding to the data set; d ij For data point γ i and γ j The distance between i For S and γ i The distance between them is less than d c The number of data points;

[0119] Step S2012: Design express A descending subscript order of q1 ≥ρ q2 ≥…≥ρ qη ; Calculate the distance as follows:

[0120]

[0121] In formula (2), Represents data points With data points The distance between When the local density is the largest, δ qi Indicates that the data set is The data point with the largest distance The distance between them; in all local densities greater than Among the data points, The data point with the smallest distance the distance between them;

[0122] Step S2013, determining the cluster center, which includes: selecting points with high local density and relatively high distance as the cluster center, assigning the remaining points to the corresponding nearest neighbor cluster with higher density, and determining the optimal mode number K = 3. The clustering result is as follows: Figure 5 As shown;

[0123] Step S202: performing noise reduction based on the VMD algorithm, which specifically includes:

[0124] The vibration signal extracted in step S1 is decomposed and denoised according to the optimal modal number K obtained in step S2013, wherein the denoising result is as follows: Figure 6 As shown, the Pearson correlation between the components after VMD decomposition and the original signal is analyzed as the basis for noise reduction and reconstruction. Figure 7 As shown;

[0125] Step S203, using the autoregressive power spectrum to determine a reasonable boundary segmentation frequency, which includes:

[0126] The AR model spectrum is used instead of the Fourier spectrum to determine the power spectrum peak and the boundary division frequency of AMD. The AR model represents the signal as follows:

[0127]

[0128] In formula (3), p is the order of the AR model; u(n) has a mean of 0 and a variance of White noise sequence; a k (k=1,2,...,p) is the corresponding p-order model parameter;

[0129] The standard equation of the autoregressive model is as follows:

[0130]

[0131] In formula (4), r x (m) is the autoregressive variable; r k (k) is the autoregressive variable sequence;

[0132] Then the autoregressive power spectrum P of the signal x(n) is AR (e jω ) can be calculated by the following formula:

[0133]

[0134] The AR model spectrum is estimated using the Burg method.

[0135] Step S204: using VGWO to determine the optimal split frequency and performing modal decomposition based on AMD, which includes:

[0136] Step S2041: Determine the first D-order damped natural frequencies of the structure ω=[ω1, ω2, ..., ω D ] T ;

[0137] Step S2042: Use VGWO to optimize the boundary segmentation frequency in AMD to obtain the optimal vibration attenuation curve, and create a population D consisting of n particles based on VGWO = (D1, D2, ..., Dn ), where D i =[d i1 ,d i2 ,...,d iD ] T , then the boundary division frequency of AMD can be expressed as:

[0138]

[0139] In formula (6), d ij is the boundary dividing frequency cutoff bandwidth of the j-th mode of the structure in the i-th particle; n is the number of particles;

[0140] Step S2043: In order to optimize the parameters using the performance index, the fitness value can be designed to be the signal x obtained by AMD decomposition and separation. j (t) and the remaining signal x k The correlation coefficient of (t) is calculated and the fitness function satisfies:

[0141]

[0142] In formula (7), ρ j is x j (t) and x k Pearson correlation coefficient of (t); and They are samples x j (t) and x k (t) mean; T is the number of samples; if x j (t) is completely separated, then x j (t) and x k The correlation coefficient of (t) is the smallest;

[0143] Step S2044: Introduce the speed component of the particle swarm optimization algorithm into the gray wolf optimization algorithm to form a variable speed gray wolf optimization algorithm. The speed and position components of the fused particle swarm algorithm are as follows:

[0144]

[0145] p i (m+1)=p i (m)+v i (m+1) (9)

[0146] In formulas (8) and (9), vi is the speed of the i-th gray wolf; p i is the current position of the i-th gray wolf; c1, c2, c3 are learning factors and satisfy c1, c2, c3∈[0,1]; ζ is the inertia factor; m is the number of iterations; X GWO1 、X GWO2 、XGWO3 Gray Wolf GWO Relative to Gray Wolf α GWO , β GWO , δ GWO The forward vector of .

[0147] Step S205, identifying modal damping, modal frequency, and modal vibration shape, specifically includes:

[0148] After decomposing the modal components of each order by AMD in step S204, the modal vibration shape is obtained by considering the structural modal response at all l sensor positions, where the k-th modal vibration shape vector of the structure is as follows:

[0149]

[0150] In formula (10), u k,j (t m ) is the kth modal response at the jth sensor position; t m is the time when the local minimum or local maximum of the modal response occurs; l is the number of installed vibration sensors; Nor(·) represents the normalized calculation; take k = 1, 2, ..., K respectively, and repeat the calculation of formula (10) to obtain the modal vibration matrix

[0151] The identification results of the modal frequency and damping ratio of the overall test bench are shown in Table 1:

[0152] Table 1 Modal frequency and damping ratio identification results of the overall test bench

[0153] Order Frequency (Hz) Damping ratio (%) 1 52.61 2.31 2 254.19 2.46 3 295.11 1.61

[0154] Step S3: For the rotor system simulation model obtained in step S2, a response surface model correction algorithm based on VGWO-DBSCAN-SVM is used to correct it, and a high-precision rotor system simulation model is constructed according to the corrected model parameters;

[0155] Specifically, in this embodiment, step S3 includes:

[0156] Step S301: determining the optimal modal number K based on the improved DBSCAN clustering algorithm;

[0157] Step S302: Fit the response surface using a support vector machine to establish an SVM response surface model between the correction variable Γ and the modal frequency. The SVM response surface model is specifically as follows: Figure 8 As shown, according to the optimal modal number obtained in step S301, the number of response surface fittings is set to K, and l sample sets of the correction variable Γ and the modal frequency are given:

[0158] {(Γ1,Θ 1,k ),(Γ2,Θ 2,k ),...,(Γ l ,Θ l,k )} (11)

[0159] In formula (11), l is the number of samples; Γ i is the variable to be corrected; Θ i,k is the sample output; k = 1, 2, ..., K is the k-th modal frequency; with the fitting accuracy ε as the error function, the optimization problem is as follows:

[0160]

[0161] In formula (12), α i and is the Lagrange coefficient; κ(Γ i ,Γ j ) is any kernel function that satisfies the Mercer condition;

[0162] Solving equation (12) yields the nonlinear regression fitting function of the kth-order modal frequency and the correction variable Γ:

[0163]

[0164] In this formula (13), b is a nonlinear term;

[0165] Step S303: using the VGWO algorithm to perform model correction, which specifically includes:

[0166] Based on the support vector machine response surface established in step S302, the measured modal frequency and the theoretical modal frequency Θ i,k The objective function is established as follows:

[0167]

[0168] In formula (14), Γ L and Γ U are the upper and lower limits of the correction variable, respectively;

[0169] More specifically, the model modification results are as follows Figure 9 The modal frequency correction results are shown in Table 2:

[0170] Table 2 Modal frequency correction results of the overall test bench

[0171]

[0172] Step S4: performing friction experiments and initializing the dynamic and static contact stiffness of the bearing rollers on the bearing system geometry simulation model to build an initial simulation model of the bearing system, wherein the dynamic and static friction coefficients of the bearings are obtained through the friction experiments;

[0173] Specifically, in this embodiment, in step S4, a friction experiment is performed to obtain the dynamic and static friction coefficient of the bearing. Figure 10 As shown, the sampling frequency is set to 200ms, the test duration is set to 10min., the load is set to 10N, the bearing is lubricated, and the average sliding friction is taken as 0.001.

[0174] Step S5: constructing a digital twin model of the rotor-bearing system based on the high-precision rotor system simulation model obtained in step S3 and the initial bearing system simulation model obtained in step S4;

[0175] Specifically, in this embodiment, step S5 includes:

[0176] In the ADAMS simulation software, the outer ring and inner ring of the bearing in the initial simulation model of the bearing system obtained in step S4 are respectively connected to the bearing seat in the high-precision rotor system simulation model obtained based on step S3 and to the rotating shaft through external connection points; at the same time, the bearing roller is set to contact with the inner and outer rings, a revolute pair is added to the rotating shaft, and a thin layer unit is added to the platform base and fixed to the ground.

[0177] Step S6: Run the rotor-bearing system test bench and the digital twin model established in step S5 under the same complex variable operating conditions, extract the output vibration signal characteristics, and perform similarity evaluation using a correlation evaluation criterion;

[0178] Specifically, in this embodiment, step S6 includes:

[0179] Step S601: performing vibration signal feature extraction, which includes:

[0180] Step S6011: For the rotor-bearing system test bench and the digital twin model, the load and speed of the rotor-bearing system are varied simultaneously. The load and speed variation range should cover the operating range of the rotor-bearing system test bench. As many operating conditions as possible are studied. Time domain vibration signals under complex variable operating conditions are extracted at the same location on the test bench and the digital twin model using the same sampling frequency.

[0181] Step S6012: Extract the obtained variable working condition vibration data respectively:

[0182] Time domain characteristics such as mean, absolute mean, effective value, average power, RMS value, peak value, peak-to-peak value, variance, standard deviation, skewness, kurtosis, peak index, waveform index, pulse index, margin index, skewness index and kurtosis index. For details of the above characteristics, see Figure 11 ;

[0183] Frequency domain features such as centroid frequency, mean square frequency, root mean square frequency, mean frequency, frequency standard deviation and frequency variance. For details of the above features, see Figure 12 ;

[0184] Wavelet energy ratio and other time-frequency domain features, the above features are specifically referred to Figure 13 ;

[0185] Step S6013: construct a feature set using the above features;

[0186] Step S602: Execute feature similarity evaluation, which specifically includes:

[0187] The time domain features, frequency domain features, and time-frequency domain features of the test bench and the digital twin model obtained in step S601 are similarly evaluated using the correlation evaluation criteria. The correlation evaluation criteria are constructed as follows:

[0188]

[0189] In formula (15), F h and l h Represent the characteristic value of the h-th sample and the corresponding time sequence respectively; and are the sample eigenvalue sequence and the time series mean, respectively; H is the number of samples; the value of the correlation evaluation index ranges from 0 to 1. The better the characteristic correlation between the test bench and the digital twin model vibration signal, the closer the value is to 1, otherwise it is closer to 0.

[0190] Step S7: For the similarity evaluation result obtained in step S6, if the similarity of the vibration signal characteristics is lower than a preset value, the dynamic and static contact stiffness of the bearing roller in the simulation model is adjusted until the similarity of the output vibration signal characteristics is higher than the preset value, and the current bearing roller contact stiffness is output. The bearing roller contact stiffness value within the full operating range of the test bench is obtained by fitting the bearing roller contact stiffness values ​​under multiple working conditions.

[0191] Specifically, in this embodiment, step S7 specifically includes:

[0192] Step S701: Obtain the dynamic and static contact stiffness of the bearing roller under a single working condition, which specifically includes:

[0193] When the similarities of all features obtained in step S602 exceed 0.7, it indicates that the vibration signal characteristics of the test bench and the digital twin model are basically consistent, and the current dynamic and static contact stiffness of the bearing roller is considered to be the bearing roller contact stiffness under the corresponding working condition;

[0194] Otherwise, the dynamic and static contact stiffness of the bearing rollers is adjusted until all characteristic similarities of the output vibration signals exceed 0.7;

[0195] Step S702: acquiring the dynamic and static contact stiffness of the bearing rollers in the full operating range, which specifically includes:

[0196] The dynamic and static contact stiffness of the bearing rollers obtained under multiple working conditions of load and speed changes are fitted to obtain the bearing contact stiffness change curve within the full operating range of the rotor-bearing test bench. The curve is as follows: Figure 14 shown.

[0197] The specific calculation results of bearing stiffness under different loads are shown in Table 3:

[0198] Table 3 Calculation results of contact stiffness between roller and inner and outer rings under different loads

[0199]

[0200] In summary, the present invention provides a method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model. For the rotor-bearing system test bench, simulation models of the rotor system and the bearing system are established respectively. For the rotor system, an algorithm based on VGWO-DBSCAN-AMD is designed to identify modal parameters and a response surface algorithm based on VGWO-DBSCAN-RF is designed to correct the simulation model, thereby obtaining a high-precision simulation model of the rotor system; for the bearing system, a friction test is carried out to obtain the friction coefficient of the bearing roller and build an initial simulation model of the bearing system; a digital twin model is obtained by assembling the simulation models of the bearing system and the rotor system. The system test bench and the digital twin model are operated under the same working conditions and the time domain vibration signals are extracted at the same time. The vibration characteristics are analyzed for correlation using the correlation evaluation criteria. The roller contact stiffness is adjusted to obtain the true value of the contact stiffness of the test bench under different working conditions. The contact stiffness under multiple working conditions is fitted to obtain the bearing roller contact stiffness value within the full operating range of the test bench.

[0201] Anything not described in detail in the present invention is well known to those skilled in the art.

[0202] The above describes in detail the preferred embodiments of the present invention. It should be understood that those skilled in the art can make numerous modifications and variations based on the concepts of the present invention without inventive effort. Therefore, any technical solutions that can be derived by those skilled in the art through logical analysis, reasoning, or limited experimentation based on the concepts of the present invention and the prior art should be within the scope of protection defined by the claims.

Claims

1. A method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model, characterized in that: The method includes: Step S1: constructing a bearing system geometry simulation model and a rotor system geometry simulation model for a rotor-bearing system test bench, and then performing a hierarchical modal test on each component in the rotor system geometry simulation model; Step S2: Based on the results of the hierarchical modal test in step S1, the modal parameters are identified using the VGWO-DBSCAN-AMD algorithm, and the modal parameter values ​​are input into the ADAMS finite element simulation model to construct a rotor system simulation model; wherein the method includes: firstly determining the optimal modal number based on the improved DBSCAN clustering algorithm, then performing noise reduction processing on the vibration signal based on the optimal modal number and the VMD algorithm, then determining the boundary splitting frequency using the autoregressive power spectrum, then determining the optimal splitting frequency using VGWO and performing modal decomposition based on AMD to obtain modal components of each order, and finally obtaining the modal vibration shape; Step S3: For the rotor system simulation model obtained in step S2, a response surface model correction algorithm based on VGWO-DBSCAN-SVM is used to correct it, and a high-precision rotor system simulation model is constructed according to the corrected model parameters; Step S4: performing friction experiments and initializing the dynamic and static contact stiffness of the bearing rollers on the bearing system geometry simulation model to build an initial simulation model of the bearing system, wherein the dynamic and static friction coefficients of the bearings are obtained through the friction experiments; Step S5: constructing a digital twin model of the rotor-bearing system based on the high-precision rotor system simulation model obtained in step S3 and the initial bearing system simulation model obtained in step S4; Step S6: Run the rotor-bearing system test bench and the digital twin model established in step S5 under the same complex variable operating conditions, extract the output vibration signal characteristics, and perform similarity evaluation using a correlation evaluation criterion; Step S7: For the similarity evaluation result obtained in step S6, if the similarity of the vibration signal characteristics is lower than a preset value, the dynamic and static contact stiffness of the bearing roller in the simulation model is adjusted until the similarity of the output vibration signal characteristics is higher than the preset value, and the current bearing roller contact stiffness is output. The bearing roller contact stiffness value within the full operating range of the test bench is obtained by fitting the bearing roller contact stiffness values ​​under multiple working conditions.

2. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 1, characterized in that: The step S1 comprises: Step S101, constructing a bearing system geometric simulation model, which includes: establishing a bearing system geometric simulation model based on ADAMS for the bearing cage, rollers, and inner and outer rings of the rotor-bearing system test bench; Step S102, constructing a rotor system geometric simulation model, which includes: constructing an ADAMS-based rotor system geometric simulation model by removing components other than the components used to construct the ADAMS-based bearing system geometric simulation model for the rotor-bearing system test bench; Step S103: performing modal tests on each component in the rotor system geometric simulation model and the rotor system itself.

3. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 1, characterized in that: In step S2, determining the optimal modal number based on the improved DBSCAN clustering algorithm specifically includes: Step S2011: Calculate the vibration response data set obtained based on the layered modal test in step S1 The local density of the data points in is as follows: In formula (1), d c represents the cutoff distance; η is the number of data sets; I S ={1,2,…,η} is the indicator set corresponding to the data set; d ij For data point γ i and γ j The distance between i For S and γ i The distance between them is less than d c The number of data points; Step S2012: Design express A descending subscript order of q1 ≥ρ q2 ≥…≥ρ qη ; Calculate the distance as follows: In formula (2), Represents data points With data points The distance between When the local density is the largest, δ qi Indicates that the data set is The data point with the largest distance The distance between them; in all local densities greater than Among the data points, The data point with the smallest distance the distance between them; Step S2013, determining the cluster center, which includes: selecting points with high local density and relatively high distance as the cluster center, assigning the remaining points to the corresponding nearest neighbor cluster with higher density, and determining the optimal mode number K=3.

4. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 3 is characterized in that: In step S2, the noise reduction process is performed on the vibration signal based on the optimal modal number and the VMD algorithm, which specifically includes: Decomposing and denoising the vibration signal extracted in step S1 according to the optimal modal number K obtained in step S2013; The Pearson correlation between the components after VMD decomposition and the original signal is analyzed as the basis for denoising and reconstruction.

5. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 4, characterized in that: In step S2, the use of the autoregressive power spectrum to determine the boundary segmentation frequency specifically includes: The AR model spectrum is used instead of the Fourier spectrum to determine the power spectrum peak and the boundary division frequency of AMD. The AR model represents the signal as follows: In formula (3), p is the order of the AR model; u(n) has a mean of 0 and a variance of White noise sequence; a k is the p-order model parameter, k=1,2,...,p; The standard equation of the autoregressive model is as follows: In formula (4), r x (m) is the autoregressive variable; r k (k) is the autoregressive variable sequence; Then the autoregressive power spectrum P of the signal x(n) is AR (e jω ) can be calculated by the following formula: The AR model spectrum is estimated using the Burg method.

6. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 5, characterized in that: In step S2, determining the optimal split frequency using VGWO and performing modal decomposition based on AMD specifically includes: Step S2041: Determine the first D-order damped natural frequencies of the structure ω=[ω1, ω2, ..., ω D ] T ; Step S2042: Use VGWO to optimize the boundary segmentation frequency in AMD to obtain the optimal vibration attenuation curve, and create a population D consisting of n particles based on VGWO = (D1, D2, ..., D n ), where D i =[d i1 ,d i2 ,...,d iD ] T , then the boundary division frequency of AMD can be expressed as: In formula (6), d ij is the boundary frequency cutoff bandwidth of the j-th mode of the structure in the i-th particle; n is the number of particles; Step S2043: In order to optimize the parameters using the performance index, the design fitness value can be obtained by decomposing and separating the signal x from AMD. j (t) and the remaining signal x k The correlation coefficient of (t) is calculated and the fitness function satisfies: In formula (7), ρ j is x j (t) and x k Pearson correlation coefficient of (t); and They are samples x j (t) and x k (t) mean; T is the number of samples; if x j (t) is completely separated, then x j (t) and x k The correlation coefficient of (t) is the smallest; Step S2044: Introduce the speed component of the particle swarm optimization algorithm into the gray wolf optimization algorithm to form a variable speed gray wolf optimization algorithm. The speed and position components of the fused particle swarm algorithm are as follows: p i (m+1)=p i (m)+v i (m+1) (9) In formulas (8) and (9), v i is the speed of the i-th gray wolf; p i is the current position of the i-th gray wolf; c1, c2, c3 are learning factors and satisfy c1, c2, c3∈[0,1]; ζ is the inertia factor; m is the number of iterations; X GWO1 、X GWO2 、X GWO3 Gray Wolf GWO Relative to Gray Wolf α GWO , β GWO , δ GWO The forward vector of The obtaining of the modal vibration shape specifically includes: After decomposing the modal components of each order by AMD in step S204, the modal vibration shape is obtained by considering the structural modal response at all l sensor positions, where the k-th modal vibration shape vector of the structure is as follows: In formula (10), u k,j (t m ) is the kth modal response at the jth sensor position; t m is the time when the local minimum or local maximum of the modal response occurs; l is the number of installed vibration sensors; Nor(·) represents the normalized calculation; take k = 1, 2, …, K respectively, and repeat the calculation of formula (10) to obtain the modal vibration matrix 7. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 6, characterized in that: The step S3 specifically includes: Step S301: determining the optimal modal number K based on the improved DBSCAN clustering algorithm; Step S302: Use support vector machine to fit the response surface to establish an SVM response surface model between the correction variable Γ and the modal frequency. The number of response surface fittings is set to K according to the optimal mode number obtained in step S301. Given l sample sets of correction variable Γ and modal frequency: {(C1,I) 1,k ),(C2,I 2,k ),…,(C l ,I l,k )} (11) In formula (11), l is the number of samples; Γ i is the variable to be corrected; Θ i,k is the sample output; k = 1, 2, …, K is the k-th modal frequency; with the fitting accuracy ε as the error function, the optimization problem is as follows: In formula (12), α i and is the Lagrange coefficient; κ(Γ i ,Γ j ) is any kernel function that satisfies the Mercer condition; Solving equation (12) yields the nonlinear regression fitting function of the kth-order modal frequency and the correction variable Γ: Step S303: using the VGWO algorithm to perform model correction, which specifically includes: Based on the support vector machine response surface established in step S302, the measured modal frequency and the theoretical modal frequency Θ i,k The objective function is established as follows: In formula (14), Γ L and Γ U are the upper and lower limits of the correction variable, respectively.

8. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 7, characterized in that: The step S5 specifically includes: By using ADAMS simulation software, the outer ring and inner ring of the bearing in the initial simulation model of the bearing system obtained in step S4 are respectively connected to the bearing seat in the high-precision rotor system simulation model obtained based on step S3 and to the rotating shaft through external connection points; at the same time, the bearing rollers are set to contact with the inner and outer rings, a revolute pair is added to the rotating shaft, and a thin layer unit is added to the platform base and fixedly connected to the ground.

9. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 8, characterized in that: The step S6 specifically includes: Step S601: performing vibration signal feature extraction, which includes: Step S6011: For the rotor-bearing system test bench and the digital twin model, the load and speed of the rotor-bearing system are simultaneously varied. The load and speed variation range should cover the operating range of the rotor-bearing system test bench. Time domain vibration signals under complex variable operating conditions are extracted at the same location on the test bench and the digital twin model using the same sampling frequency. Step S6012: extracting the following features from the acquired variable working condition vibration data: time domain features, frequency domain features, and time-frequency domain features; Step S6013: construct a feature set using the above features; Step S602: Execute feature similarity evaluation, which specifically includes: The time domain features, frequency domain features, and time-frequency domain features of the test bench and the digital twin model obtained in step S601 are similarly evaluated using the correlation evaluation criteria. The correlation evaluation criteria are constructed as follows: In formula (15), F h and l h Represent the characteristic value of the h-th sample and the corresponding time sequence respectively; and are the sample eigenvalue sequence and the time series mean, respectively; H is the number of samples; the value of the correlation evaluation index ranges from 0 to 1. The better the characteristic correlation between the test bench and the digital twin model vibration signal, the closer the value is to 1, otherwise it is closer to 0.

10. The method for measuring the dynamic and static contact stiffness of bearing rollers based on a digital twin model according to claim 9, characterized in that: The step S7 specifically includes: Step S701: Obtain the dynamic and static contact stiffness of the bearing roller under a single working condition, which specifically includes: When the similarities of all features obtained in step S602 exceed 0.7, it indicates that the vibration signal characteristics of the test bench and the digital twin model are basically consistent, and the current dynamic and static contact stiffness of the bearing roller is considered to be the bearing roller contact stiffness under the corresponding working condition; Otherwise, the dynamic and static contact stiffness of the bearing rollers is adjusted until all characteristic similarities of the output vibration signals exceed 0.7; Step S702: acquiring the dynamic and static contact stiffness of the bearing rollers in the full operating range, which specifically includes: The dynamic and static contact stiffness of the bearing rollers obtained under multiple working conditions with varying loads and speeds are fitted to obtain the bearing contact stiffness variation curve within the full operating range of the rotor-bearing test bench.

Images