Rubber vibration isolator reliability evaluation method based on digital twinning

Abstract

This invention discloses a reliability assessment method for rubber vibration isolators based on digital twins, belonging to the field of reliability assessment technology. The method includes the following steps: conducting accelerated aging tests on the rubber material of the isolator, collecting degradation data of the rubber material at multiple time points, fitting the evolution of the rubber elastic modulus with aging time and temperature using a four-parameter Arrhenius degradation model to obtain the prior distribution of the four-dimensional parameters; constructing a high-fidelity finite element digital twin model of the isolator structure, establishing a deterministic mapping from the material-level elastic modulus to the structural-level dynamic stiffness, extracting the simulated dynamic stiffness at a preset calibration frequency point and the measured dynamic stiffness at the corresponding frequency point in the structural-level high-temperature vibration coupling test; and using a higher-order moment method to approximate high-dimensional weighted integration through sparse deterministic nodes to construct a posterior degradation model, and then assessing the reliability of the isolator under arbitrary operating conditions. This invention can provide an effective assessment of the state perception and remaining life of high-reliability rubber vibration isolation products.

Description

Technical Field

[0001] This invention relates to the field of vibration isolator structural reliability assessment, and specifically to a method for assessing the reliability of rubber vibration isolators based on digital twins. Background Technology

[0002] Vibration isolators in aero-engine systems are critical components connecting the engine and the airframe structure, and their performance directly affects the vibration transmission characteristics and structural lifespan of the aircraft. Typical rubber vibration isolators consist of a metal-rubber composite structure, combining rubber materials with a metal frame through a vulcanization process. They utilize the energy storage characteristics of rubber elasticity to absorb and isolate vibration energy, and their reliability directly impacts the operational stability of the power system, passenger comfort, and the service life of critical components.

[0003] Under actual service conditions, rubber vibration isolators are subjected to high-temperature environments and broadband random vibration excitation for extended periods, leading to continuous thermo-oxidative aging and cyclic stress fatigue of the rubber material. Consequently, the dynamic stiffness of the isolator gradually deviates from the design value over service time. Dynamic stiffness, as a core indicator characterizing the force-displacement relationship of the isolator under forced vibration, directly affects the vibration transmission characteristics of the isolation system. Its deviation from the design value will cause frequency band drift and uncontrolled vibration energy transmission path, thereby interfering with the measurement accuracy of airborne precision instruments, accelerating fatigue failure of key components such as engine bearings, affecting the stable operation of spacecraft in orbit, and even triggering low-frequency resonance cracking of welded structures.

[0004] Therefore, it is necessary to conduct dynamic stiffness degradation modeling and reliability assessment of rubber vibration isolators under complex working conditions in order to avoid economic losses or irreversible catastrophic accidents. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention proposes a high-order moment-Bayes update-maximum entropy reliability assessment framework based on the concept of digital twins. A four-parameter Arrhenius degradation model is used to describe the evolution of the rubber elastic modulus with aging time and temperature, and parameter priors are identified based on accelerated aging data of material-grade dumbbell specimens. The high-order moment method generates 81 deterministic key sampling points and their corresponding weights in a 4-dimensional parameter space. A high-fidelity finite element model is used to physically map the elastic modulus to structural-level dynamic stiffness at each sampling point. The Bayesian likelihood is efficiently evaluated using sparse quadrature, and the structural-level dynamic stiffness data obtained from multi-stress reliability enhancement tests is used to sequentially update the parameter posterior. On the updated rubber degradation model, the uncertainty carried by the parameter posterior is further efficiently mapped to the dynamic stiffness response under arbitrary temperature-time conditions using the high-order moment method, obtaining its first three original moments. The probability density function of the dynamic stiffness is reconstructed using the maximum entropy principle, and the reliability of the vibration isolator is evaluated accordingly. The proposed framework can efficiently integrate material-level prior and structural-level measured information under the condition of scarce structural-level data, balancing computational efficiency and statistical accuracy.

[0006] Specifically, the present invention provides a reliability assessment method for rubber vibration isolators based on digital twins, which includes the following steps: S1. Conduct accelerated aging tests on the rubber material of the vibration isolator and collect data on the degradation of the rubber material at multiple time points during the aging process; S2. The evolution of the rubber elastic modulus with aging time and temperature is fitted using a four-parameter Arrhenius degradation model to obtain a rubber elastic modulus degradation model. Combined with the rubber material degradation data collected at multiple time points, a four-dimensional parameter multivariate Gaussian prior is obtained. S3. Construct a high-fidelity finite element digital twin model of the vibration isolator structure, establish a deterministic mapping from material-level elastic modulus to structural-level dynamic stiffness, and extract the simulated dynamic stiffness at the preset calibration frequency point. S4. Conduct high-temperature vibration coupling tests on the structure of the vibration isolator and extract the measured dynamic stiffness at the same preset calibration frequency point in the test. S5. Based on a four-dimensional parameter multivariate Gaussian prior, a sequential update strategy of alternating measured dynamic stiffness at two temperatures is implemented. A high-order moment method is used to approximate high-dimensional weighted integrals through sparse deterministic nodes. Material-level prior data are corrected using structural-level measured data to construct a posterior degradation model that closely resembles the actual service state.

[0007] S6. Based on the a priori degradation model, the reliability of the vibration isolator under arbitrary operating conditions is evaluated using the high-order moment maximum entropy method.

[0008] Preferably, S2 includes: S21. Establish a four-parameter Arrhenius degradation model to describe the evolution of the elastic modulus with aging time and temperature, perform a logarithmic transformation on the parameters, and convert E... a / R is used as the overall parameter to obtain the four-dimensional parameter vector θ; S22. Using the collected rubber material degradation data from multiple time points as input, and assuming that the measurement noise is independently and identically distributed according to a normal distribution, a problem of minimizing the sum of squared residuals is constructed and solved using the LM algorithm to obtain parameter point estimates. residual vector The Jacobian matrix J; S23. Based on the asymptotic normality of nonlinear least squares, For mean vector, For the covariance matrix, construct a four-dimensional parameter multivariate Gaussian prior.

[0009] Preferably, the rubber elastic modulus degradation model in S2 is: (1); Where E0 is the unaged elastic modulus, Z is the pre-exponential factor, and E α Let α be the apparent activation energy, R be the gas constant, α be the time power law exponent, and T be the absolute temperature.

[0010] Preferably, S3 includes: S31. Based on the actual geometry and service conditions of the vibration isolator, a three-dimensional solid model consisting of rubber elements, aluminum alloy frame and counterweight mass block is established in the three-dimensional software and imported into the finite element preprocessing module. S32. Set the element size according to the geometric characteristics and stress gradient of each component, divide the hexahedral mesh, and set the rubber-metal vulcanization interface and the contact between the counterweight and the upper end face as the bonded contact to obtain the finite element mesh model. S33. An isotropic linear elastic constitutive model is adopted for the rubber element, and its elastic modulus is provided by the degenerate model for each working condition. A constant parameter linear elastic constitutive model is adopted for the metal skeleton to obtain a finite element model that imparts material properties. S34. Set boundary conditions according to the actual installation method of the test bench; S35. Using harmonic response analysis, a vertical simple harmonic excitation force is applied to the upper surface, and the steady-state displacement amplitude at the preset calibration frequency is extracted to obtain the simulated dynamic stiffness.

[0011] Preferably, the deterministic mapping from material-level elastic modulus to structural-level dynamic stiffness in S3 is as follows: ; in, (·) denotes the finite element operator defined by the harmonic response process, and E(t,T) is the rubber elastic modulus degradation model.

[0012] Preferably, S4 includes: S41. Conduct high-temperature and random vibration coupling tests on the vibration isolator; S42. Arrange piezoelectric accelerometers on the upper surface of the counterweight and at the lower foundation of the vibration isolator to obtain the power spectral density and calculate the acceleration transmissibility. S43. Convert the transmissivity into dynamic stiffness and extract the measured dynamic stiffness at the preset calibration frequency point.

[0013] Preferably, S5 includes: S51. Using the obtained four-dimensional parameters as the initial prior, establish a Gaussian measurement noise model for structural dynamic stiffness observation and construct the likelihood function. The posterior distribution of the parameters is given by Bayes' theorem. S52. The structural dynamic stiffness observation data under dual temperatures are sorted alternately according to the aging time to form a dual-temperature alternating observation sequence, and the posterior parameters are sequentially recursively extrapolated point by point according to the observation sequence. S53. In each sequential update step, the current prior is mapped to the four-dimensional standard independent normal space through Cholesky decomposition. The one-dimensional Gaussian quadrature node and weight are constructed using the higher-order moment MQ method, and then expanded by tensor product to obtain sparse deterministic four-dimensional quadrature node and corresponding weight. S54. Transform the four-dimensional quadrature nodes back to the parameter space to obtain prior parameter samples. Use the degenerate model and finite element operator to obtain the predicted dynamic stiffness values ​​of each node. Combine the measured dynamic stiffness to calculate the likelihood values ​​of each node and normalize them to obtain the posterior weights. S55. The posterior mean vector and posterior covariance matrix are obtained by weighted summation of each node with posterior weights and used as the prior for the next observation point. After traversing all observation points, the final posterior distribution of parameters is obtained, and a posterior degradation model that closely approximates the actual service state is constructed.

[0014] Preferably, the posterior degradation model in S5 is: ; ; ; ; Where μ0 is the mean vector and Σ0 is the covariance matrix. Let K be the measured dynamic stiffness at the k-th observation point. The simulated dynamic stiffness at the k-th observation point when the parameter is θ. To measure noise, To measure the noise variance, Let K be the likelihood function, and K be the total number of structural observation points.

[0015] Preferably, S6 includes: S61. Propagate parameter uncertainties to elastic modulus E via a rubber elastic modulus degradation model; S62. Perform finite element mapping by MQ node arrangement and extract the dynamic stiffness at the preset calibration frequency. S63. The first three orders of the dynamic stiffness at the origin are obtained analytically by summing the nodal responses using the Gaussian quadrature formula. S64. Reconstruct the dynamic stiffness probability density from the first three original moments of the dynamic stiffness based on the maximum entropy. S65. Reliability is calculated based on the probability that the dynamic stiffness probability density is greater than the threshold.

[0016] Preferably, the method further includes: S7. Repeat the high-order moment maximum entropy method evaluation process on the mesh where temperature and aging time are coupled to obtain the dynamic stiffness degradation zone and reliability degradation curve at each temperature.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: The high-order moment-Bayesian-maximum entropy framework provided by this invention integrates material-level priors from accelerated aging of rubber samples with structural-level data obtained from multi-stress reliability enhancement tests to assess the reliability of vibration isolators. It achieves cross-level information complementarity through Bayesian theory, compensates for insufficient structural test samples using low-cost material aging tests, and reduces testing costs. It avoids model bias caused by artificially predetermined probability distributions by combining the maximum entropy method, and achieves multi-level uncertainty propagation analysis using the high-order moment method. This scheme takes into account both the intrinsic aging mechanism of rubber materials and the macroscopic structural failure characteristics of vibration isolators, and can obtain more objective and accurate life and reliability assessment results under limited structural test data conditions. It is suitable for accelerated aging reliability analysis and durability evaluation of rubber vibration isolators.

[0018] This invention propagates multi-source uncertainties via a high-fidelity finite element physics engine using a higher-order moment method, obtaining higher-order statistical moments of dynamic stiffness without relying on traditional Monte Carlo sampling methods. This offers four advantages: First, it significantly reduces computational costs. Traditional Monte Carlo methods require a large number of random samples and iterative solutions to dynamic stiffness. This invention establishes a mapping from random parameters to higher-order moments of dynamic stiffness through a non-sampling approach, greatly reducing the number of repeated solutions to frequency-domain dynamic stiffness, making it particularly suitable for uncertainty analysis of multi-frequency scanning and multi-degree-of-freedom structures. Second, it improves solution stability. Monte Carlo higher-order moment estimation is significantly affected by sample size and prone to random fluctuations; the statistical moments obtained by this invention are deterministic, free from sampling disturbances, and effectively avoid distortions in higher-order moments such as skewness and kurtosis under finite samples. Third, it expands applicable scenarios. It can better handle problems such as drastic changes in dynamic stiffness near resonance and non-normal distribution of random responses. The obtained higher-order statistical moments can be directly combined with methods such as maximum entropy to reconstruct the probability density function, facilitating the quantification of structural vibration uncertainties and reliability analysis. Fourth, it facilitates iterative engineering applications. The analytical framework is computationally controllable, easy to embed into structural optimization and uncertain parameter inversion cyclic calculations, and has good engineering application value.

[0019] 2. This invention proposes a sparse Bayesian likelihood weighting method based on aPC high-order moment nodes, replacing the high-dimensional sampling of traditional MCMC with 81 4-dimensional deterministic nodes, and adopting a dual-temperature staggered sequential update strategy to synchronously inject information from the two temperatures into the parameter space. The experimental results show that after only 16 steps of staggered updates, the standard deviation of each parameter is significantly reduced, and the prediction accuracy of the validation set is significantly improved compared with the prior.

[0020] 3. On the response prediction side, this invention obtains the first three moments of dynamic stiffness using 5-node one-dimensional higher-order moments and high-fidelity finite element simulation, and reconstructs its probability density function using the maximum entropy principle. In a single working condition, only 5 finite element calls are needed to complete the reliability assessment, which significantly reduces the computational threshold of probability analysis and can quantitatively characterize the temperature-sensitive reliability degradation law of rubber vibration isolators at different temperatures.

[0021] 4. This invention further reveals the physical fact that the actual service aging rate is faster than the material-level Arrhenius extrapolation through the difference analysis of the a priori / posterior elastic modulus surfaces, and verifies the core value of the method of this invention in "continuously correcting the virtual model with real data" through digital twin. Attached Figure Description

[0022] Figure 1 This is a flowchart of the reliability assessment method for rubber vibration isolators based on digital twins according to the present invention; Figure 2 This is a graph showing the degradation of the elastic modulus of rubber at 100°C and 120°C as a function of aging time, according to an embodiment of the present invention. Figure 3This is a three-dimensional schematic diagram of the vibration isolator structure and mass block according to an embodiment of the present invention; Figure 4 This is a finite element mesh diagram of the vibration isolator structure and mass block according to an embodiment of the present invention; Figure 5 This is an acceleration power spectral density diagram of the lower surface of a rubber elastic support according to an embodiment of the present invention; Figure 6 This is a comparison diagram of the prior and posterior four-parameter distributions of an embodiment of the present invention; Figure 7 This is a comparison diagram of the measured dynamic stiffness, prior predicted dynamic stiffness, and posterior predicted dynamic stiffness at a verification point at 100°C, according to an embodiment of the present invention. Figure 8 This is a comparison chart of the verification errors of prior predicted dynamic stiffness and posterior predicted dynamic stiffness at 100°C in an embodiment of the present invention. Figure 9 This is a comparison diagram of the measured dynamic stiffness, prior predicted dynamic stiffness, and posterior predicted dynamic stiffness at a verification point at 120°C, according to an embodiment of the present invention. Figure 10 This is a comparison chart of the verification errors of prior predicted dynamic stiffness and posterior predicted dynamic stiffness at 120°C in an embodiment of the present invention. Figure 11 This is a dynamic stiffness degradation curve of a ±3σ band at 100°C according to an embodiment of the present invention. Figure 12 This is a dynamic stiffness degradation curve of a ±3σ band at 110°C according to an embodiment of the present invention. Figure 13 This is a dynamic stiffness degradation curve of a ±3σ band at 120°C according to an embodiment of the present invention. Detailed Implementation

[0023] Hereinafter, embodiments of the present invention will be described with reference to the accompanying drawings.

[0024] As attached Figure 1 As shown, as an example of an implementation of the present invention, a rubber vibration isolator for a certain type of aircraft engine is used as the evaluation object. The isolator has external dimensions of 145×145×90mm and is made of natural rubber elastomer and 7075 aluminum alloy skeleton by vulcanization bonding. The entire device weighs approximately 6.77kg. The following is a complete demonstration of the entire process from material-level aging data, through degradation modeling, finite element mapping, Bayesian updating to reliability assessment, using actual data. The changes in data between each step and the final result can be clearly seen.

[0025] This invention provides a reliability assessment method for rubber vibration isolators based on digital twins, comprising the following steps: S1. Conduct accelerated aging tests on the rubber material of the vibration isolator and collect data on the degradation of the rubber material at multiple time points during the aging process.

[0026] To obtain the degradation law of the mechanical properties of rubber materials with aging time and temperature, and to provide basic data for subsequent degradation model fitting, accelerated aging tests of rubber materials were carried out.

[0027] As an embodiment of the present invention, natural rubber from the same batch as the vibration isolator was used to prepare dumbbell-shaped standard specimens, and accelerated thermo-oxidative aging tests were conducted under a constant compression ratio of 0.3 (replicating the pre-compression state in service). Two temperature levels, 100℃ and 120℃, were set, with a total aging time of 18 hours. Samples were taken at 10 time points: 0, 1, 2, 3, 4, 6, 8, 10, 14, and 18 hours. The stress-strain curves of each specimen were collected using a universal testing machine, and the elastic modulus E was calculated. The average of parallel specimens was taken to obtain 20 elastic modulus observation data at the two temperatures, which served as the initial raw input data for this method, as shown in Table 1.

[0028] Table 1. Accelerated Aging Test Data of Material-Grade Elastic Modulus (Compression Ratio 0.3)

[0029] As shown in Table 1, the elastic modulus of rubber decreases monotonically with aging time, and the decay is significantly faster at 120℃ than at 100℃, demonstrating the temperature sensitivity of thermo-oxidative aging to the mechanical properties of rubber. The slight scattering in the data reflects measurement noise, which forms the basis for subsequent least-squares fitting and statistical estimation of the noise variance. The degradation curves of rubber elastic modulus with aging time (at 100℃ and 120℃) are shown below. Figure 2 As shown.

[0030] S2. A four-parameter Arrhenius degradation model is used to fit the evolution of the rubber elastic modulus with aging time and temperature. Combined with rubber material degradation data collected at multiple time points, a four-dimensional multivariate Gaussian prior is obtained, specifically including: S21. Establish a four-parameter Arrhenius degradation model to describe the evolution of the elastic modulus E(t,T) with aging time and temperature, perform a logarithmic transformation on the parameters, and convert E... a / R is used as the overall parameter to obtain the four-dimensional parameter vector θ.

[0031] The chemical aging of rubber in a thermo-oxidative environment is mainly dominated by two processes: molecular chain breakage and re-crosslinking. The evolution of its macroscopic properties is commonly described using the Arrhenius model. This invention employs a four-parameter Arrhenius degradation model to describe the evolution of the rubber's elastic modulus E(t,T) with aging time t and temperature T: (1); Where E0 is the unaged elastic modulus in MPa; Z is the pre-exponential factor in h. -α E α The apparent activation energy is expressed in J·mol⁻¹. -1 R = 8.314 J·mol -1 ·K -1 Let α be the gas constant, α be the time power law exponent, and T be the absolute temperature in K.

[0032] Introduce logarithmic transformations ln E0 and ln Z, and then... α / R is used as the overall parameter, resulting in a four-dimensional parameter vector θ: θ = [ ln E0, ln Z , E α / R , α ] (2).

[0033] S22. Using the collected rubber material degradation data from multiple time points as input, and assuming that the measurement noise is independently and identically distributed according to a normal distribution, a problem of minimizing the sum of squared residuals is constructed and solved using the Levenberg-Marquardt algorithm to obtain parameter point estimates. residual vector With Jacobian matrix J.

[0034] In order to obtain the parameters of the degradation model from the material-level test data, the present invention uses the least squares method for fitting.

[0035] Assume that the measurement noises are independent and all follow the same normal distribution. Then the maximum likelihood estimate of the parameters is equivalent to minimizing the sum of squared residuals: ; in, Least squares (maximum likelihood) estimation of the four-dimensional parameter vector θ; N m The total number of material-level elastic modulus observation data (N in this embodiment) m =20); Let i be the measured elastic modulus; When the parameter is θ, the aging time of the degradation model at the i-th observation point is given by equation (1). absolute temperature The given value is the predicted value of the elastic modulus; argmin represents the parameter θ that minimizes the sum of squared residuals.

[0036] This invention uses the Levenberg-Marquardt algorithm (LM algorithm, damped least squares method) to solve equation (3). Substituting the observed data, the parameter point estimates are obtained. and the corresponding residual vector and Jacobian matrix .

[0037] S23. Based on the asymptotic normality of nonlinear least squares, For mean vector, Given the covariance matrix, construct a four-dimensional parameter multivariate Gaussian prior p0(θ)=N(μ0,Σ0).

[0038] Due to measurement noise variance Knowing this beforehand, this invention performs unbiased estimation based on the mean square of residuals: ; in, For measuring the unbiased estimate of the noise variance; For residual vectors 2-norm square (sum of squared residuals); N m -4 represents the degrees of freedom (number of observations minus the number of parameters to be estimated, which is 4).

[0039] Furthermore, based on the asymptotic normality of nonlinear least squares, parameter estimation... The true value is obtained by following the mean and the covariance is... The approximate multivariate Gaussian distribution. The multivariate Gaussian prior with this as a parameter in this invention:

[0040] .

[0041] Among them, the mean This refers to the parameter estimates obtained from the fitting, specifically the covariance. The Jacobian matrix reflects the uncertainty in parameter identification caused by measurement noise. This characterizes the sensitivity of each parameter to the model output.

[0042] As an embodiment of the present invention, the 20 observation data in Table 1 are substituted into (Equations 3-5), assuming that the measurement noise is independently and identically distributed in N(0, σ). E 2 The Levenberg-Marquardt algorithm is used to solve the problem of minimizing the sum of squared residuals, and the parameter point estimates are obtained. residual vector The Jacobian matrix J.

[0043] Noise variance estimated from residual mean square unbiased Then, by applying the asymptotic normality of nonlinear least squares, we obtain the multivariate Gaussian prior of the parameters p0(θ) = N(μ0, Σ0), where .

[0044] The marginal prior distributions of the parameters obtained from the solution are as follows: ln E0~ N(1.6074, (2.09×10 -2 ) 2 ) ln Z ~ N(11.335, 3.26 2 ) E α / R ~ N(5546.7, (1.30×10 3 ) 2 ) ; α ~ N(0.5837, (9.76×10) -2 ) 2 ) Where, N(μ, σ) 2 ) represents a mean of μ and a variance of σ. 2 It follows a normal distribution.

[0045] This prior corresponds to the unaged elastic modulus E0 = e 1.6074 ≈ 4.99 MPa, apparent activation energy E α = (E α / R)·R≈ 5546.7×8.314 ≈ 46.1 kJ·mol -1 .

[0046] This is the initial estimate of the degradation parameters that can be obtained solely from material-level data, which serves as the prior input for subsequent Bayesian updates—thus transforming the data from 20 discrete elastic modulus points into a probability distribution of four-dimensional parameters.

[0047] S3. Construct a high-fidelity finite element digital twin model of the vibration isolator structure, establish a deterministic mapping from material-level elastic modulus to structural-level dynamic stiffness, and extract the simulation dynamic stiffness at preset calibration frequency points, specifically including: S31. Based on the actual geometry and service conditions of the vibration isolator, a three-dimensional solid model consisting of rubber elements, an aluminum alloy frame and a counterweight mass block is established in the three-dimensional software and imported into the finite element preprocessing module.

[0048] The multi-scale mapping from material-level elastic modulus to structural-level dynamic stiffness is the physical engine of the digital twin "virtual body". Based on the actual geometry and service conditions of rubber elastic supports, this invention establishes a high-fidelity three-dimensional solid finite element model in a finite element platform as the physical engine of the digital twin, realizing a deterministic mapping from rubber elastic modulus to structural-level dynamic stiffness.

[0049] The rubber elastic support of the research object has external dimensions of 145 × 145 × 90 mm and consists of two parts: a rubber component and a 7075 aluminum alloy bracket, which are bonded together using a vulcanization process. To realistically reproduce the load-bearing state of the vibration isolator under service conditions, a geometric model was created on top of the elastic support using a 1:1 scale steel counterweight actually used in the test setup, and this counterweight was fixed to the upper surface of the elastic support. The overall model, including the rubber component, aluminum alloy bracket, and steel counterweight, was imported into the finite element preprocessing module after being modeled in SolidWorks.

[0050] S32. Set the element size according to the geometric characteristics and stress gradient of each component, divide the hexahedral mesh, and set the rubber-metal vulcanization interface and the contact between the counterweight and the upper end face as the bonded contact to obtain the finite element mesh model.

[0051] Based on the actual geometry and service conditions of the vibration isolator, a three-dimensional solid model consisting of a rubber element, a 7075 aluminum alloy frame, and a steel mass block was created in the three-dimensional parametric mechanical CAD software SolidWorks (the mass block was modeled at a 1:1 scale of the test bench and fixed to the upper end face of the vibration isolator with bolts). A three-dimensional schematic diagram of the vibration isolator structure and the mass block is shown below. Figure 3 As shown, import the data into the finite element preprocessing module. The rubber element is meshed with a 5mm mesh, the aluminum alloy frame with a 10mm mesh, and the mass block with a 50mm mesh. The finite element mesh diagrams of the vibration isolator structure and the mass block are shown below. Figure 4 As shown, the entire model comprises 26,568 solid elements and 43,360 nodes. The rubber-metal vulcanization interface and the contact between the counterweight and the upper surface are both set as bonded contacts. The model is meshed using hexahedral elements, with element sizes set according to the geometric characteristics and stress gradients of each component: the rubber component, as the main deformation and energy dissipation carrier, has an element size of 5mm; the aluminum alloy support has a gentler stress gradient, so the element size is 10mm; and the counterweight has an element size of 50mm to ensure stress transfer accuracy in the assembly area. The mesh yielded 26,568 solid elements and 43,360 nodes. For interface treatment, the rubber-metal vulcanization interface uses bonded contact to simulate the strong bonding state under vulcanization adhesion; the counterweight and the upper surface of the elastic support are mechanically connected by nuts. Because the preload applied by the nuts prevents separation between the bottom surface of the counterweight and the upper surface of the support under the investigated small-amplitude vibration excitation, this interface also uses bonded contact. The resulting mesh is shown below. Figure 4 As shown: S33. An isotropic linear elastic constitutive model is adopted for the rubber element, and its elastic modulus is provided by the degenerate model for each working condition. A constant parameter linear elastic constitutive model is adopted for the metal skeleton to obtain a finite element model that imparts material properties.

[0052] Within the investigated service temperature range and under small-amplitude vibration excitation, the working strain level of the rubber component is relatively low, falling within the near-linear elastic range. Therefore, this invention does not introduce hyperelastic constitutive models such as Mooney-Rivlin, but instead directly characterizes the mechanical properties of rubber using the elastic modulus E, employing an isotropic linear elastic constitutive model. The rubber elastic modulus E at each temperature and aging time is given by the S2 rubber elastic modulus degradation model and updated into the finite element model for each operating condition. The metal support material is 7075 aluminum alloy, which also adopts a linear elastic constitutive model, and the material parameters are considered constant within the investigated service temperature range.

[0053] As an embodiment of the present invention, the rubber is in the near-linear elastic region under the condition of small-amplitude vibration superimposed with pre-compression, so it is directly characterized by the elastic modulus E (provided by the degradation model of Equation (1) for each working condition at each temperature and aging time). An isotropic linear elastic constitutive model is adopted in combination with geometric nonlinearity (large deformation) to capture the geometric stiffening caused by pre-compression. The initial parameters of the two materials are shown in Table 2.

[0054] Table 2 Initial Material Parameters

[0055] S34. Set boundary conditions according to the actual installation method of the test bench; Boundary conditions were set according to the actual installation method of the vibration table: a fixed support constraint was applied to the surface of the bolt holes at the bottom of the support, equivalent to simulating the rigid connection between the bottom of the launcher and the foundation; the counterweight mass was fixed to the upper surface of the support by bolt fasteners, which was realized in the finite element model as a combination of concentrated mass and preload, with a top preload of 6000 N to simulate the static preload effect of the engine weight on the elastic support in actual service. The preload was applied in a static step before solving the dynamic stiffness, and the resulting stress-deformation field served as the initial configuration for subsequent harmonic response analysis to ensure that the geometric stiffening effect under preload conditions was correctly transmitted.

[0056] S35. Using harmonic response analysis, a vertical simple harmonic excitation force is applied to the upper surface, and the steady-state displacement amplitude |x(f0)| at the preset calibration frequency is extracted to obtain the simulated dynamic stiffness.

[0057] This invention uses dynamic stiffness K as the output quantity to reflect the overall force-displacement relationship of an elastic support structure under forced vibration. Considering that rubber elastic supports mainly bear low-frequency vibration excitation in actual service, and that 6.3 Hz is the center frequency of the 1 / 3 octave band specified in ISO 266, a commonly used standard frequency point in low-frequency vibration testing and dynamic characteristic evaluation of rubber components; previous frequency sweep results also show that the displacement response amplitude of this elastic support is most significant at 6.3 Hz, which is the most representative frequency of mechanical response during service. Therefore, this invention uses 6.3 Hz as the evaluation frequency point for dynamic stiffness, defining dynamic stiffness as the ratio of the excitation force amplitude to the steady-state displacement response amplitude at this frequency point, ensuring consistency between simulation and experimental results.

[0058] Using the harmonic response analysis module of finite element software, a simple harmonic excitation force with an amplitude of F0 = 1 N was applied vertically to the upper surface of the elastic support structure. Fixed support constraints were applied at the bottom of the elastic support. The frequency sweep range covered the main response frequency band of the structure, and the excitation point on the upper surface was extracted. steady-state displacement amplitude at The dynamic stiffness at that frequency point is obtained by the ratio of the force amplitude F0 to the displacement amplitude under harmonic excitation. ; in, The simulated dynamic stiffness at the preset calibration frequency; The amplitude of the simple harmonic excitation force is F0 (in this embodiment, F0 = 1N); |x(f0)| is the amplitude of the steady-state displacement response at the excitation point at frequency f0; f0 is the preset calibration frequency (in this embodiment, f0 = 6.3 Hz).

[0059] Combining the above geometry, mesh, material, boundary conditions, and solution process, the established finite element model constitutes a system from the rubber elastic modulus E to the structural dynamic stiffness. Deterministic mapping: ; in, (·) represents the finite element operator defined by the harmonic response process of equation (5), and E(t,T) is the rubber elastic modulus degradation model given by equation (1).

[0060] Therefore, the finite element model consists of the rubber elastic modulus E and the structural dynamic stiffness D. f Deterministic mapping D f = F(E).

[0061] After calibration, the mapping is approximately linear near the operating point, with a sensitivity of approximately 2.94 kN·mm. -1 / MPa: For example, the initial dynamic stiffness D is obtained by mapping the unaged elastic modulus E0 = 4.93MPa (posterior value) through finite element method. f,0=38.35 kN·mm -1 .

[0062] This finite element model, or digital twin physical model, serves as the unified computational basis for "predicting dynamic stiffness from parameters" in subsequent Bayesian updates and "predicting dynamic stiffness from modulus" in reliability assessments.

[0063] S4. Conduct high-temperature vibration coupling tests on the structure of the vibration isolator, and extract the measured dynamic stiffness at the same preset calibration frequency point in the test, specifically including: S41. Conduct high-temperature and random vibration coupling tests on the vibration isolator.

[0064] A high-temperature vibration coupling test was conducted on a certain type of rubber elastic support structure. The overall dimensions are 145 mm × 145 mm × 90 mm, and the mass is about 6.77 kg. It is made of a 7075 aluminum alloy shell and a natural rubber elastomer vulcanized and bonded together.

[0065] According to standards such as GJB 150.3A-2009 High Temperature Test Standard and GJB 150.16A-2009 Vibration Test Standard, random vibration coupling tests at 100 ℃ and 120 ℃ were conducted. The tests were carried out in a walk-in three-dimensional integrated test chamber, with fixed support constraints applied to the bottom of the elastic support and a counterweight fixed at the top to simulate the actual service load.

[0066] S42. Arrange piezoelectric accelerometers on the upper surface of the counterweight and at the lower foundation of the vibration isolator to obtain the power spectral density and calculate the acceleration transmissibility T(f).

[0067] This invention, based on GB 50463-2019 "Engineering Vibration Isolation Design Standard", defines the failure criterion as the dynamic stiffness degradation rate and uses dynamic stiffness K as the intrinsic mechanical parameter. According to national standard GB / T 22159.3-2008, this invention uses an indirect method to extract the dynamic stiffness of the rubber elastic support. During the test, a counterweight of known mass m is fixed above the elastic support, and piezoelectric accelerometers are placed on the upper surface of the counterweight and at the lower foundation of the support. The acceleration time history signals at the output and input ends are obtained through a charge amplification and data acquisition system. and The acquired signal is frequency-domain transformed to obtain the acceleration power spectral density at both ends. and The actual measurement results of the two are as follows Figure 5 As shown.

[0068] Depend on and The square root of the ratio is defined as the acceleration transmissibility: .

[0069] S43. Convert the transmissivity into dynamic stiffness and extract the measured dynamic stiffness at the preset calibration frequency point.

[0070] Within the effective frequency band where the counterweight acts as a rigid mass and its excitation frequency is significantly lower than its internal resonant frequency, a frequency domain relationship can be established between the inertial force on the counterweight and the relative displacement between the two ends of the elastic support. The dynamic stiffness can then be calculated using the transmissibility. ; extract The dynamic stiffness is denoted as K, with units of kN / mm. This extraction method is physically consistent with the dynamic stiffness defined in finite element simulation, facilitating the modeling of structural-level experimental data and material-level priors. The measured dynamic stiffness K extracted from equation (31) is used as the measured dynamic stiffness in equation (10). .

[0071] Due to limitations in test duration and cost, only a limited number of observation points were obtained in the early degradation stage. The measured data of structural-level dynamic stiffness degradation are shown in Table 3.

[0072] Table 3. Measured Degradation Data of Structural Dynamic Stiffness

[0073] Note: The last four time points marked with "*" (t=8.5, 10, 13, 14h) are independent validation sets; the measured sequences of the first eight time points (t=0, 1.5, 2.5, 3, 4, 5, 6, 7.5h) are used for Bayesian updates. Thus, a total of 12 cumulative aging time points and 12 dynamic stiffness observations at each of the two temperatures are used as "physical observation" inputs to correct the material-level priors.

[0074] S5. Execute a dual-temperature alternating Bayesian sequential update strategy, employing the high-order moment MQ method to approximate a high-dimensional weighted integral through sparse deterministic nodes. Use structural-level measured data to correct material-level prior data, constructing a posterior degradation model that closely approximates the actual service state. Specifically, this includes: S51. Using the four-dimensional parameters obtained in S2 as the initial prior, establish a Gaussian measurement noise model for structural dynamic stiffness observation and construct the likelihood function. The posterior distribution of the parameters is given by Bayes' theorem.

[0075] Based on material-level identification results, the degradation model parameter vector The prior is taken as a 4-dimensional joint normal distribution: ; in, Let S2 represent a multivariate normal distribution with mean vector μ0 and covariance matrix Σ0, where μ0 and Σ0 are given by S2.

[0076] Suppose that the measured values ​​of the structural dynamic stiffness of the vibration isolator were obtained at k observation points. .remember This represents the predicted dynamic stiffness value given by the finite element digital model under the corresponding working conditions.

[0077] Assuming that the measurement noises are independent and identically distributed according to a zero-mean normal distribution, then: ; in, The measured dynamic stiffness at the k-th observation point; =ℱ(E(t_k,T_k;θ)) is the predicted dynamic stiffness of the finite element digital model under the working condition (t_k,T_k) when the parameter is θ; For measuring noise; To measure the noise variance, it is estimated from the root mean square residuals of the structural-level measured sequence and the finite element baseline sequence under the prior mean. The corresponding likelihood function is: ; in, is the likelihood function; K is the total number of structural observation points.

[0078] By Bayes' theorem, the posterior distribution of the parameters satisfies: ; in, p0(θ) is the posterior distribution of the parameter; p0(θ) is the prior of equation (9).

[0079] because It is a solid element finite element operator that considers the elasticity and geometric nonlinearity of rubber. The cost of a single call is not negligible. Directly estimating the posterior by sampling methods such as MCMC on equation (12) in the 4-dimensional parameter space usually requires 10 3 -10 4 The sheer scale of finite element method calls is prohibitive in engineering practice.

[0080] Therefore, this invention employs the higher-order moment method to perform sparse deterministic node approximation of the high-dimensional weighted integral of the molecule in equation (12), requiring only The Gaussian quadrature estimate of the posterior mean vector and the posterior covariance matrix can be given in one step by a single finite element call. Its node construction, tensor product extension and single-step Bayesian weighting formula are shown in Equations (14)-(29).

[0081] S52. The structural dynamic stiffness observation data under dual temperatures are sorted alternately according to the aging time to form a dual-temperature alternating observation sequence, and the posterior parameters are sequentially recursively extrapolated point by point according to this order.

[0082] To fully utilize the complementary information of structural dynamic stiffness at two temperature levels, 100 ℃ and 120 ℃, this invention employs a dual-temperature staggered sequential update strategy: updating the dynamic stiffness of the structure at both temperatures... The observation points are arranged alternately in chronological order, resulting in a length of... observation sequence The two adjacent steps correspond to observations at 100 ℃ and 120 ℃ at the same time, respectively. Under this sequence, the posterior parameter is recursively derived using the following formula: ; in, To absorb the posterior distribution after the kth observation (as the prior of the (k+1)th step); , The aging time and temperature are for the k-th observation; This represents the number of observation points used for updating at a single temperature. For material-level priors, This represents the likelihood for a single observation point.

[0083] Each step uses the posterior of the previous step as the prior, accumulating observation information point by point; the staggered order of the two temperatures ensures that their contributions are synchronously injected into the parameter space, avoiding posterior drift caused by "using one temperature first and then using another temperature".

[0084] S53. In each sequential update step, the current prior is mapped to a four-dimensional standard independent normal space through Cholesky decomposition. One-dimensional Gaussian quadrature nodes and weights are constructed using the higher-order moment method, and sparse deterministic four-dimensional quadrature nodes and corresponding weights are obtained through tensor product expansion.

[0085] The posterior recursion of equation (13) and the subsequent reliability assessment based on the degradation model both involve the weighted integration of high-fidelity finite element operators in a high-dimensional parameter space; directly using sampling methods such as MCMC results in excessive computation.

[0086] The basic principles of the higher-order moment method are first introduced below. Then, its specific application in the two stages of Bayesian posterior update and posterior response moment calculation in this invention is given. Finally, the dynamic stiffness probability density function is reconstructed in combination with the maximum entropy principle to support reliability assessment.

[0087] The core idea of ​​the Moment Quadrature (MQ) method is to construct an orthogonal polynomial basis of the input probability density and corresponding Gaussian quadrature node-weight pairs based solely on the finite-order primitive moments of the input random variables. This allows the propagation process from the input distribution to the output statistical moments to be analytically weighted from the model responses at a small number of key sampling points, without the need for large-scale random sampling. The specific construction steps are as follows: Let the input random variable ξ take values ​​in the interval Ω, with probability density p(ξ), and its first 2k raw moments be: ; Construct using the first 2k orders of the original moments Hankel matrix of order: ; Cholesky decomposition exists in positive time. ,in Let be an upper triangular matrix. By Mysovskih's theorem, we can construct a theorem for the weight function. A number of orthogonal polynomials: ; This family of orthogonal polynomials satisfies three recurrence relations: ; The recursive coefficient is: ; By compactly rewriting the three recursive terms of equation (17) into matrix form, we can obtain the Jacobi matrix. Its diagonal element is The second diagonal element is By applying the Golub-Welsch theorem, eigenvalue decomposition of the Jacobi matrix yields k Gaussian quadrature nodes and their corresponding weights: ; in, for The first component of the corresponding unit eigenvector, Let be the j-th eigenvalue of the Jacobi matrix J. The quadrature formula derived from this is: ; For any polynomial of degree not exceeding 2k-1 Strictly true, therefore only The result can be obtained from the model evaluation. The higher-order moments are the moments of each order, which is why the higher-order moment method obtains higher-order statistics with very few model calls.

[0088] In each step of the posterior recursion of equation (13) above, what is actually needed is not the posterior density itself, but several posterior statistics, all of which can be written in the form of posterior expectations. Directly calculating such expectations requires high-dimensional integration of the posterior over a 4-dimensional parameter space, but the posterior contains unnormalized constants, making direct sampling difficult. Note that by Bayes' theorem... The posterior expectation can be strictly rewritten as the ratio of the two prior expectations: ; in, [·]、 [·] denotes the expected value of the posterior and prior distributions, respectively; g(θ) is the function corresponding to any statistic to be determined.

[0089] Specifically, let the prior at step k be a multivariate normal distribution. By using Cholesky decomposition, the parameter space is mapped to a 4-dimensional standard independent normal space. : ; Where ξ~N(0,I4) is a four-dimensional standard independent normal random vector; For the current prior covariance The Cholesky factorization of the lower triangular factors satisfies .

[0090] This transformation not only achieves decorrelation, but also The coupling relationship carried is completely preserved in In the spatial node layout.

[0091] In standardized space In the middle, for each component, take a one-dimensional standard normal MQ. Each node has 4 dimensions, and the weights are formed by the tensor product of the 1-dimensional nodes. ; in, , (d=1,...,4;j_d=1,..., () represents the nodes and weights of the one-dimensional standard normal higher-order moments in the direction of the d-th component; This represents the number of nodes per dimension.

[0092] This invention takes ,common Each 4D MQ node can rigorously and accurately integrate a 5th-order polynomial in each component direction, and its engineering precision is sufficient to cover the requirements for low-order moment propagation under a 4D Gaussian prior.

[0093] S54. Transform the four-dimensional quadrature nodes back into the parameter space to obtain prior parameter samples. Use the degenerate model and finite element operator to obtain the predicted dynamic stiffness values ​​of each node. Combine the measured dynamic stiffness to calculate the likelihood values ​​of each node and normalize them to obtain the posterior weights.

[0094] The 81 standardized nodes of equation (23) are transformed back to the parameter space by equation (22) to obtain 81 parameter samples in the prior space. These samples are then substituted into equations (1)-(2) to calculate the corresponding elastic modulus. The 81 predicted dynamic stiffness values ​​are obtained by the finite element operator. ; Where j is the four-dimensional node number, j=1,...,81; Let j be the j-th standardized four-dimensional node; This is the j-th parameter sample obtained by inverse transformation back to the parameter space using equation (22); Take parameters The time-dependent degradation model (1) under operating conditions The elastic modulus is given below; The finite element operator defined by equation (8); Let be the predicted dynamic stiffness value of the j-th node.

[0095] Based on the Gaussian likelihood of equation (11) at a single observation point, the likelihood values ​​at each node are: ; From the Gaussian quadrature formula of equation (20), any a priori function The expectation can be approximated by a weighted sum of the nodes: ; Substitute equation (26) into the numerator of equation (21) respectively (take... ) and denominator (take By using the denominator as a normalization constant, the normalized posterior weights on the nodes can be defined: .

[0096] S55. The posterior mean vector and posterior covariance matrix are obtained by weighted summation of each node with posterior weights and used as the prior for the next observation point. After traversing all observation points, the final posterior distribution of parameters is obtained, and a posterior degradation model that closely approximates the actual service state is constructed.

[0097] Any posterior expectation can be derived from the posterior weights. The weighted summation at the prior nodes directly yields: ; Specifically, in equation (28), take respectively and You can get the first one. The posterior mean vector and the posterior second-order raw moment are used; the posterior covariance matrix is ​​obtained by subtracting the outer product of the posterior mean from the latter. ; Higher-order moment method sequential update (dual temperature alternation, prior → posterior): Starting with the multivariate Gaussian prior p0(θ)=N(μ0,Σ0) of S2, the first 8 time points of S4, totaling 16 structural-level dynamic stiffness observation points at both temperatures, are absorbed into the parameter space point by point. To fully utilize the complementary information of the two temperatures, a dual-temperature alternating sequential update strategy is adopted: at the same time, the observations at 100℃ are absorbed first, followed by the observations at 120℃, avoiding posterior drift caused by exhausting the information at a certain temperature first.

[0098] The likelihood function takes the single-point Gaussian form: L(θ|D r ) = ∏ (1 / (√(2π)σ m )) · exp( -(D r,k - D f,k (θ)) 2 / (2σ m 2 ) ) ; Among them, D f,k (θ)=F(E(t k ,T k ;θ)) represents the finite element method's prediction of dynamic stiffness under the corresponding working condition; the standard deviation of the measurement noise σ m The value is taken as a fixed constant and is estimated by the root mean square of the residuals between the structural-level measured sequence and the finite element baseline sequence under the prior mean parameters.

[0099] By Bayes' theorem p(θ|D) r ) ∝ p0(θ)·L(θ|D r ).

[0100] Since F(·) is a finite element operator for solid elements, the cost of a single evaluation is high. Directly using MCMC to sample in the four-dimensional parameter space requires 10... 3 -10 4 Sub-finite element evaluation is impractical in engineering. Therefore, a higher-order moment MQ method is used to approximate a high-dimensional weighted integral with sparse deterministic nodes: P1, at each step, the current prior N(μ) is... k-1 , Σ k-1 (Decomposed by Cholesky) Mapped to a four-dimensional standard independent normal space .

[0101] P2, Each component takes n GP =3 one-dimensional standard normal MQ nodes, tensor product forms N Q =3 4 =81 four-dimensional nodes and their weights.

[0102] P3. Transform the 81 standardized nodes back into the parameter space to obtain 81 prior parameter samples, and substitute them into equation (1) to obtain the corresponding elastic modulus.

[0103] P4. 81 dynamic stiffness predictions D obtained from finite element operators. f,i =F(E(t k ,T k ;θ i )).

[0104] P5. Calculate the likelihood value L for each node. i Define normalized posterior weights, and then directly obtain the posterior mean from the weighted sum of the nodes. With posterior covariance And serve as a prerequisite for the next step.

[0105] Posterior mean and covariance:

[0106] ; .

[0107] Each step requires only 81 finite element evaluations, yielding μ k Σ k That is, it serves as a priori for the next step.

[0108] The evolution trajectory of the four parameters in the 16-step sequential update process (i.e. how the parameter data is gradually corrected with observation) is shown in Table 4.

[0109] Table 4. Two-Temperature Alternating Sequential Bayesian Update Trajectory Table

[0110] After 16 updates, the posterior parameters converged to ln E0 = 1.596, ln Z = 11.109, and E α / R=5357.8, α=0.502, corresponding to E0=4.93 MPa, E α =44.54 kJ·mol -1 .

[0111] like Figure 6 As shown, compared to the prior, the posterior distributions of all four parameters significantly shrink and undergo a center shift: α shifts left from the prior mean of 0.584 to 0.502 (the structural-level data substantially corrects the long-term shape of the degradation curve); E α / R decreased from 5546.70 to 5357.83, and the standard deviation decreased from 1.30 × 10 3 Shrinking to 1.26×10 3ln E0 shifted slightly from 1.607 to 1.596; ln Z remained basically near its mean, mainly showing a further compression of the distribution width.

[0112] At this point, the material-level prior data has been corrected by structural-level observations into a posterior degradation model that more closely resembles actual service conditions.

[0113] Furthermore, the posterior prediction is validated to verify the improvement in data accuracy: The four verification time points (t=8.5, 10, 13, 14h) that were not involved in the update were substituted into the prior model and the posterior model respectively, and compared with the measured dynamic stiffness at the structural level. The comparison of prediction errors is shown in Table 5.

[0114] Table 5. Comparison of Prediction Errors on the Validation Set (Unit: kN·mm) -1 )

[0115] At 100℃, the prior model's predictions are systematically overestimated (mean error 0.3563 kN·mm). -1 After Bayesian updating, the mean error of the posterior model was reduced to 0.0983 kN·mm. -1 The accuracy was improved by 72.4%, such as... Figure 7 , Figure 8 As shown in the figure, the gray dashed line represents the residual between the predicted dynamic stiffness and the measured dynamic stiffness; at 120℃, the average error is 0.6405 kN·mm. -1 Reduced to 0.3400 kN·mm -1 The accuracy was improved by 46.9%, such as Figure 9 , Figure 10 As shown in the figure, the gray dashed line represents the residual between the predicted dynamic stiffness and the measured dynamic stiffness.

[0116] This indicates that structural-level dynamic stiffness measurements carry information not available in material-level priors (the true modulus degradation trajectory of structural components under multiaxial non-uniform stress states, and the cyclic stress-thermal-oxygen coupling aging effect). The higher-order moment method Bayesian update framework incorporates these into the parameters identified by material-level specimens, achieving effective cross-scale correction from material-level priors to structural-level posteriors.

[0117] S6. Based on the a priori degradation model, the reliability of the vibration isolator under arbitrary (t,T) working conditions is evaluated using the high-order moment maximum entropy method.

[0118] The following example demonstrates the complete data flow from parameter posteriority to reliability, using the service conditions T=120℃ and t=50h as an example.

[0119] S61. The parameter uncertainty is propagated to the elastic modulus through the rubber elastic modulus degradation model.

[0120] Because the degenerate model has weak local nonlinearity of the parameters, after first-order (delta method) approximation, E under the target condition approximately follows a normal distribution with a mean μ. E Standard deviation σ E The posterior mean and covariance of the parameters are analyzed. Substituting T = 120℃ = 393.15K and t = 50h, the posterior parameters are given as follows: E ~ N( 2.7757 , 0.154 2 ) MPa (mean 2.7757 MPa, standard deviation 0.154 MPa).

[0121] S62. Perform finite element mapping through MQ node arrangement to extract the dynamic stiffness at the preset calibration frequency.

[0122] MQ node deployment → Finite element mapping: in E ~ N(μ E , σ E 2 Five key sampling points ζ are arranged using the higher-order moment method. j and corresponding weight w j (i.e., 5 Gauss-Hermitian nodes), each E j Substitute each point into the finite element model of S3 to solve for the harmonic response, and extract the dynamic stiffness D at 6.3Hz. f,j =F(E j ).

[0123] The results are shown in Table 6. It can be seen that each modulus sample is mapped to a dynamic stiffness value by the physical model.

[0124] Table 6. Key nodes of MQ and their finite element dynamic stiffness response (T=120℃, t=50h)

[0125] S63. The first three orders of the dynamic stiffness at the origin are obtained analytically by summing the nodal responses using the Gaussian quadrature formula. Dynamic stiffness moment calculation: The first three dynamic stiffness moments are obtained analytically by weighted summation of the nodal responses using the Gaussian quadrature formula: m1 = Σw j ·D f,j =32.005, m2=Σw j ·D f,j 2 =1024.525, m3=Σw j ·D f,j 3 =32803.070, from which the average dynamic stiffness is calculated to be 32.005 kN·mm. -1 The standard deviation is approximately 0.453 kN·mm. -1The skewness is approximately 0.29.

[0126] The entire response moment calculation requires only 5 finite element evaluations, which is far less than the order of magnitude required for the Monte Carlo method to achieve the same engineering accuracy.

[0127] S64. The dynamic stiffness probability density is reconstructed from the first three original moments of the dynamic stiffness based on the maximum entropy.

[0128] The reliability cannot be directly given by moments; the dynamic stiffness probability density f(D) needs to be reconstructed from the first three moments. f ).

[0129] The maximum entropy principle, under the constraint of only knowing a finite number of moments, selects the distribution with the maximum information entropy (i.e., without introducing additional artificial assumptions) as the minimum partial estimator. Its analytical form is the exponential family f(x) = exp(Σλ). n x n ), Lagrange multiplier λ n Determined by moment constraints.

[0130] In this example, the first three moments (N=3) are taken. The third moment is introduced on the basis of the Gaussian approximation to capture the skewness and more accurately reflect the tail behavior of dynamic stiffness (which is especially important for failure probability estimation).

[0131] S65. Reliability is calculated based on the probability that the dynamic stiffness probability density is greater than the threshold.

[0132] The failure threshold is taken as a certain proportion D of the initial dynamic stiffness. f,th =η·D f,0 η=0.8.

[0133] The finite element baseline of the posterior degradation model gives the unaged dynamic stiffness D at t=0. f,0 =38.35 kN·mm -1 The corresponding failure threshold D f,th =0.8 × 38.35 = 30.68 kN·mm -1 .

[0134] Reliability is defined as the probability that the dynamic stiffness is greater than a threshold: R(t,T) = P{ D f >D f,th} = ∫ f(D f ) dD f (Integral lower limit D) f,th (to +∞); The probability density for maximum entropy reconstruction is at a threshold of 30.68 kN·mm. -1 Based on the above integration, the reliability of this operating condition is R = 0.9983.

[0135] That is, after 50 hours of thermo-oxidative aging, the average dynamic stiffness decreased from the initial 38.35 kN·mm. -1 Reduced to approximately 32.0 kN·mm -1 The reliability is 0.9983, indicating that the vibration isolator is still in a reliable working state under this operating condition.

[0136] At this point, the input "parameter posterior" data is transformed into the desired result of a single reliability value through "modulus distribution → dynamic stiffness of 5 nodes → third moment → maximum entropy probability density → threshold integration".

[0137] S7. Repeat the high-order moment maximum entropy method evaluation process on the mesh where temperature and aging time are coupled to obtain the dynamic stiffness degradation zone and reliability degradation curve at each temperature.

[0138] To reveal the reliability degradation law under the coupled effect of temperature and aging time, the complete high-order moment-maximum entropy evaluation process was repeated on a grid with T∈{100,110,120}℃ and t∈[1,300]h, and the dynamic stiffness degradation zone (mean ±3σ) and reliability degradation curve at each temperature were obtained.

[0139] Regarding dynamic stiffness, the dynamic stiffness at all three temperatures decreases monotonically with aging time, and the ±3σ band gradually widens over time (reflecting the cumulative propagation of posterior parameter uncertainty to dynamic stiffness response under long-term extrapolation).

[0140] Temperature has a significant impact on the degradation rate: degradation is slowest at 100℃, with the average value still slightly above the failure threshold after 300h; degradation accelerates significantly at 110℃, with the average value reaching the threshold at about 173h; degradation is most severe at 120℃, with the average value falling below the threshold at about 85h.

[0141] Regarding reliability, the dynamic stiffness distribution reconstructed by the maximum entropy at each moment is integrated above the failure threshold to obtain the reliability degradation curves at three temperatures. The key results are shown in Table 7.

[0142] The time for reliability to first drop to 0.9 shortens significantly with increasing temperature: 222.5 h at 100℃, 120.5 h at 110℃, and 64.3 h at 120℃. For example... Figure 11 , Figure 12 , Figure 13 As shown, all three curves exhibit a pattern of "remaining basically unchanged in the early stage and rapidly decreasing in the later stage." The higher the temperature, the shorter the high-reliability plateau and the steeper the decay slope. The time for the reliability to first drop to 0.9 shortens quasi-exponentially with increasing temperature, which is consistent with the Arrhenius acceleration law.

[0143] Table 7 Key results of reliability degradation at different temperatures

[0144] In summary, this embodiment fully demonstrates the entire data evolution process of the proposed method using actual data: Starting with aging data of 20 elastic moduli at the material level → Gaussian prior of a four-parameter degradation model is obtained through fitting → Finite element digital model is constructed to establish the mapping from E to D_f (sensitivity approximately 2.94 kN·mm). -1 ( / MPa) → Sixteen dynamic stiffness observation points at the structural level are introduced and sequentially updated using the 81-node higher-order moment method Bayesian method to obtain the posterior parameters of the shrinkage offset (ln E0=1.596, ln Z=11.109, E a (R=5357.8, α=0.502) → The prediction accuracy on the validation set is improved by 72.4% and 46.9% respectively. → The dynamic stiffness probability density is obtained by reconstruction with the 5-node high-order moment method and maximum entropy under any working condition, and the reliability is obtained by integration (e.g., R=0.9983 for 120℃ / 50h working condition). → Finally, the reliability degradation curve and reliability lifetime of temperature-time coupling are given (the time for reliability to drop to 0.9 is 222.5 / 120.5 / 64.3h at 100 / 110 / 120℃ respectively).

[0145] In terms of computational cost, only 5 finite element evaluations are required for each (t,T) condition in the reliability assessment phase, and only 81 finite element evaluations are required for each step in the update phase. This is far lower than the number of finite element simulations required for the Monte Carlo method to achieve the same accuracy in four-dimensional parameter space, demonstrating the engineering deployability of this method.

[0146] The present invention and its embodiments have been described above. This description is not restrictive, and the accompanying drawings are only one embodiment of the present invention. The actual method steps and data are not limited thereto. In conclusion, if those skilled in the art are inspired by this description and design similar method steps and embodiments without departing from the spirit of the present invention, such designs should fall within the protection scope of the present invention.

Claims

1. A reliability assessment method for rubber vibration isolators based on digital twins, characterized in that, It includes the following steps: S1. Conduct accelerated aging tests on the rubber material of the vibration isolator and collect data on the degradation of the rubber material at multiple time points during the aging process; S2. The evolution of the rubber elastic modulus with aging time and temperature is fitted using a four-parameter Arrhenius degradation model to obtain a rubber elastic modulus degradation model. Combined with the rubber material degradation data collected at multiple time points, a four-dimensional parameter multivariate Gaussian prior is obtained. S3. Construct a high-fidelity finite element digital twin model of the vibration isolator structure, establish a deterministic mapping from material-level elastic modulus to structural-level dynamic stiffness, and extract the simulated dynamic stiffness at the preset calibration frequency point. S4. Conduct high-temperature vibration coupling tests on the structure of the vibration isolator and extract the measured dynamic stiffness at the same preset calibration frequency point in the test. S5. Based on the four-dimensional parameter multivariate Gaussian prior, a dual-temperature measured dynamic stiffness alternating sequential update strategy is implemented. The high-order moment method is used to approximate the high-dimensional weighted integral through sparse deterministic nodes. The material-level prior data is corrected using structural-level measured data to construct a posterior degradation model that closely resembles the actual service state. S6. Based on the a priori degradation model, the reliability of the vibration isolator under arbitrary operating conditions is evaluated using the high-order moment maximum entropy method.

2. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 1, characterized in that, S2 include: S21. Establish a four-parameter Arrhenius degradation model to describe the evolution of the elastic modulus with aging time and temperature, perform a logarithmic transformation on the parameters, and convert E... a / R is used as the overall parameter to obtain the four-dimensional parameter vector θ; S22. Using the collected rubber material degradation data from multiple time points as input, and assuming that the measurement noise is independently and identically distributed according to a normal distribution, a problem of minimizing the sum of squared residuals is constructed and solved using the LM algorithm to obtain parameter point estimates. residual vector The Jacobian matrix J; S23. Based on the asymptotic normality of nonlinear least squares, For mean vector, For the covariance matrix, construct a four-dimensional parameter multivariate Gaussian prior.

3. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 2, characterized in that, The degradation model for the elastic modulus of rubber in S2 is as follows: （1）； Where E0 is the unaged elastic modulus, Z is the pre-exponential factor, and E α Let α be the apparent activation energy, R be the gas constant, α be the time power law exponent, and T be the absolute temperature.

4. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 1, characterized in that, S3 include: S31. Based on the actual geometry and service conditions of the vibration isolator, a three-dimensional solid model consisting of rubber elements, aluminum alloy frame and counterweight mass block is established in the three-dimensional software and imported into the finite element preprocessing module. S32. Set the element size according to the geometric characteristics and stress gradient of each component, divide the hexahedral mesh, and set the rubber-metal vulcanization interface and the contact between the counterweight and the upper end face as the bonded contact to obtain the finite element mesh model. S33. An isotropic linear elastic constitutive model is adopted for the rubber element, and its elastic modulus is provided by the degenerate model for each working condition. A constant parameter linear elastic constitutive model is adopted for the metal skeleton to obtain a finite element model that imparts material properties. S34. Set boundary conditions according to the actual installation method of the test bench; S35. Using harmonic response analysis, a vertical simple harmonic excitation force is applied to the upper surface, and the steady-state displacement amplitude at the preset calibration frequency is extracted to obtain the simulated dynamic stiffness.

5. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 4, characterized in that, The deterministic mapping from material-level elastic modulus to structural-level dynamic stiffness in S3 is as follows: ； in, (·) denotes the finite element operator defined by the harmonic response process, and E(t,T) is the rubber elastic modulus degradation model.

6. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 1, characterized in that, S4 include: S41. Conduct high-temperature and random vibration coupling tests on the vibration isolator; S42. Arrange piezoelectric accelerometers on the upper surface of the counterweight and at the lower foundation of the vibration isolator to obtain the power spectral density and calculate the acceleration transmissibility. S43. Convert the transmissivity into dynamic stiffness and extract the measured dynamic stiffness at the preset calibration frequency point.

7. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 1, characterized in that, S5 include: S51. Using the obtained four-dimensional parameters as the initial prior, establish a Gaussian measurement noise model for structural dynamic stiffness observation and construct the likelihood function. The posterior distribution of the parameters is given by Bayes' theorem. S52. The structural dynamic stiffness observation data under dual temperatures are sorted alternately according to the aging time to form a dual-temperature alternating observation sequence, and the posterior parameters are sequentially recursively extrapolated point by point according to the observation sequence. S53. In each sequential update step, the current prior is mapped to the four-dimensional standard independent normal space through Cholesky decomposition. The one-dimensional Gaussian quadrature node and weight are constructed using the higher-order moment MQ method, and then expanded by tensor product to obtain sparse deterministic four-dimensional quadrature node and corresponding weight. S54. Transform the four-dimensional quadrature nodes back to the parameter space to obtain prior parameter samples. Use the degenerate model and finite element operator to obtain the predicted dynamic stiffness values ​​of each node. Combine the measured dynamic stiffness to calculate the likelihood values ​​of each node and normalize them to obtain the posterior weights. S55. The posterior mean vector and posterior covariance matrix are obtained by weighted summation of each node with posterior weights and used as the prior for the next observation point. After traversing all observation points, the final posterior distribution of parameters is obtained, and a posterior degradation model that closely approximates the actual service state is constructed.

8. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 7, characterized in that, The posterior degeneracy model in S5 is: ； ； ； ； Where μ0 is the mean vector and Σ0 is the covariance matrix. Let K be the measured dynamic stiffness at the k-th observation point. The simulated dynamic stiffness at the k-th observation point when the parameter is θ. To measure noise, To measure the noise variance, Let K be the likelihood function, and K be the total number of structural observation points.

9. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 1, characterized in that, S6 include: S61. Propagate parameter uncertainties to elastic modulus E via a rubber elastic modulus degradation model; S62. Perform finite element mapping by MQ node arrangement and extract the dynamic stiffness at the preset calibration frequency. S63. The first three orders of the dynamic stiffness at the origin are obtained analytically by summing the nodal responses using the Gaussian quadrature formula. S64. Reconstruct the dynamic stiffness probability density from the first three original moments of the dynamic stiffness based on the maximum entropy. S65. Reliability is calculated based on the probability that the dynamic stiffness probability density is greater than the threshold.

10. The reliability assessment method for rubber vibration isolators based on digital twins according to claim 1, characterized in that, The method also includes: S7. Repeat the high-order moment maximum entropy method evaluation process on the mesh where temperature and aging time are coupled to obtain the dynamic stiffness degradation zone and reliability degradation curve at each temperature.

Images