Method for calculating dynamic stiffness of bearing seat of new energy decelerator shell

Abstract

The application belongs to the field of mechanical engineering, and particularly relates to a method for calculating dynamic stiffness of a new energy reducer shell bearing seat. The method aims to solve the problem of low dynamic stiffness prediction accuracy caused by ignoring the manufacturing residual stress and multi-physical field coupling in traditional simulation. By constructing a refined finite element model containing high-density hybrid mesh and material nonlinear characteristics, integrating the neural network prediction-based die casting residual stress field, coupling multi-source dynamic loads such as gear alternating torque, electromagnetic pulsating force and road excitation, performing nonlinear transient dynamics solving, and obtaining time-domain response data of bearing seat nodes, the application also establishes a process-structure-performance database and a fast prediction model, supports lightweight design closed-loop optimization and digital twin online monitoring, improves simulation accuracy and engineering application efficiency, and is suitable for dynamic performance analysis and life management of new energy vehicle reducers.

Description

Technical Field

[0001] This application belongs to the field of mechanical engineering, and specifically relates to a method for calculating the dynamic stiffness of the bearing seat of a new energy reducer housing. Background Technology

[0002] With the continuous advancement of new energy vehicle technology, the electric drive system, as the core power unit, directly impacts the overall vehicle's operating efficiency and reliability through its structural performance. The reducer, a key component of the electric drive system, plays a crucial role in torque transmission and speed regulation, while its housing structure must possess sufficient static strength and dynamic stability under lightweight design conditions. Among these components, the bearing housing, as the core supporting the drive shaft system, has a decisive influence on the system's modal distribution, vibration transmission path, and fatigue life. Especially under high-frequency excitation conditions, stiffness performance becomes a key indicator for evaluating the structure's dynamic response capability.

[0003] In particular, traditional engineering methods for stiffness analysis of key components of the reducer housing often employ static finite element simulation combined with empirical formula corrections. These methods typically simplify the bearing housing region to an equivalent support boundary, neglecting the nonlinear contact relationships and material anisotropy in the actual load path. This results in calculations that fail to accurately reflect the dynamic mechanical behavior under real-world conditions. Furthermore, in multiphysics coupled environments, thermal stress deformation caused by temperature gradients and alternating loads generated by electromagnetic excitation act on the housing structure, further exacerbating the time-varying and spatially non-uniform nature of the stiffness characteristics.

[0004] Current technologies generally rely on bench tests to obtain dynamic stiffness data. However, the testing process is limited by sensor placement, excitation methods, and the accuracy of boundary condition simulations, often resulting in large measurement deviations and poor repeatability. Furthermore, traditional simulation models lack the ability to explicitly model the impact of manufacturing processes (such as die-casting defects and residual stress), failing to effectively correlate microstructural features with macroscopic dynamic performance. Moreover, existing analysis workflows often treat modal analysis and harmonic response analysis separately, failing to establish an integrated computational framework from natural frequency identification to complex stiffness parameter extraction, leading to insufficient characterization of phase lag and energy dissipation mechanisms near the resonance region. Therefore, there is an urgent need for a method for calculating the dynamic stiffness of bearing seats in new energy vehicle reducers that can comprehensively consider structural nonlinearity, time-varying loads, and process disturbance effects. Summary of the Invention

[0005] The purpose of this invention is to provide a method for calculating the dynamic stiffness of the bearing seat of a new energy reducer housing, which can effectively solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows: a method for calculating the dynamic stiffness of the bearing seat of a new energy reducer housing, comprising the following specific steps: Step (1) Constructing a refined finite element model: Based on the three-dimensional geometric model of the new energy reducer housing, a high-density tetrahedral and hexahedral mixed mesh is divided, and local mesh refinement is implemented in the bearing seat area. The minimum element size is controlled within 0.5 mm, and a mesh refinement layer is set in stress concentration areas such as bolt connection holes, stiffener transition areas, and casting fillets. At the same time, anisotropic elastic modulus tensor and nonlinear constitutive relation are introduced in the material property definition. To characterize the microstructure directionality and plastic yielding behavior of die-cast aluminum alloys; Step (2) Integrate manufacturing process disturbance field: Obtain historical data stream of die-casting process parameters in actual production, establish a residual stress prediction neural network model, input variables include mold temperature, pouring speed, holding time and cooling rate, output is the three-dimensional residual stress distribution field inside the shell, use the prediction result as a prestress field superimposed on the finite element model, and combine it with X-ray diffraction measurement data to calibrate the model, with the error controlled within 8%; Step (3) Couple multi-physics dynamic load: Construct the operating conditions of electric drive system The comprehensive excitation model is applied simultaneously with the periodic alternating torque caused by gear meshing, the high-frequency axial ripple load generated by the electromagnetic pulsation of the motor, and the random vibration displacement boundary conditions transmitted from the road surface excitation to the suspension point during vehicle operation. The temperature field is imported through the thermo-mechanical coupling analysis module, and the working temperature rise range is set to 25℃ to 120℃. The thermal expansion coefficient of the material changes linearly with temperature. Step (4) performs nonlinear transient dynamic solution: using an explicit integral algorithm combined with damping correction technology, the full-time domain response calculation is completed under time step adaptive control, and the sampling frequency is not less than 20 kHz. Record the displacement response sequence and corresponding reaction force time history curve of each node on the bearing housing mounting surface in the radial, axial and tangential three degrees of freedom; Step (5) Extract complex dynamic stiffness parameters: Perform synchronous Fourier transform on the obtained displacement and reaction force signals to obtain the transfer function matrix in the frequency domain, and obtain the modal parameters of each order based on the least squares complex frequency domain method, including natural frequency, damping ratio and mode participation factor, and further calculate the complex stiffness K*=K′+iK″, where the real part K′ represents the energy storage stiffness and the imaginary part K″ represents the energy dissipation characteristics. Finally, generate the frequency-dynamic stiffness spectrum and Nyquist plot for performance evaluation.

[0007] Preferably, in step (1), the mesh generation adopts an adaptive meshing strategy based on geometric curvature and gradient change rate to ensure that the local refinement mechanism is automatically triggered in the region with curvature greater than 0.02 mm, and the size ratio of adjacent units does not exceed 1.2 times, so as to avoid numerical ill-conditioning caused by mesh distortion; in addition, the Chaboche viscoplastic model is introduced in the material nonlinear modeling, which includes three sets of back stress tensors and nonlinear hardening parameters to describe the Bauschinger effect and ratchet strain accumulation process under cyclic loading.

[0008] Preferably, in step (2), the residual stress prediction neural network adopts a hybrid architecture of convolutional long short-term memory with a depth of 6, wherein the convolutional layer is used to extract spatial correlation features, the LSTM layer captures the time series dependence of process parameters, the training set contains no less than 50,000 sets of actual production data samples, and the coefficient of determination R² after cross-validation is no less than 0.93; the residual stress field output by the model is mapped to the finite element solid mesh node in voxelized data format, and the field quantity jump phenomenon is eliminated by smooth interpolation algorithm.

[0009] Preferably, in step (3), the electromagnetic pulsating load is modeled based on the harmonic components of the motor stator winding current, mainly considering the 6th, 12th, and 18th spatial harmonic components. The amplitude is obtained by inversion using electromagnetic simulation software and converted into an equivalent simple harmonic force acting on the central axis of the bearing housing. The phase difference is dynamically adjusted based on the feedback information from the rotor position sensor. The road excitation generates a random signal based on the C-level road power spectral density function, which is then transmitted back to the reducer suspension support through the vehicle multibody dynamics model to form the displacement boundary input.

[0010] Preferably, in step (4), the explicit integration algorithm uses an improved version of the central difference method, and introduces the Hilber-Hughes-Taylor α parameter to adjust the high-frequency damping. The value of α ranges from -0.05 to -0.1, so as to suppress non-physical high-frequency oscillations without affecting the accuracy of low-frequency response. At the same time, an energy conservation monitoring module is set up. When the total fluctuation of the system's kinetic energy and internal energy exceeds 3% of the initial energy, the time step is automatically reduced and the local calculation is restarted.

[0011] Preferably, in step (5), the original signal is subjected to Hanning window weighting before Fourier transform to prevent spectral leakage and the frequency resolution is controlled within 1 Hz; the modal order range is limited to the 1st to the 15th order during the least squares complex frequency domain fitting process, and the weight function is assigned according to the signal-to-noise ratio distribution to prioritize the fitting accuracy of the resonance peak region; the complex stiffness result is normalized according to the international standard ISO 10846 and the output unit is kilonewtons per micrometer.

[0012] Preferably, it also includes: establishing a process-structure-performance related database to store the material parameters, process conditions, mesh configuration and final dynamic stiffness curve used in each calculation, forming a structured dataset; training a regression-type support vector machine model based on the dataset to quickly predict the dynamic stiffness frequency response characteristics of the new design scheme, with a prediction response time of less than 3 seconds and a relative error of less than 7%.

[0013] Preferably, it also includes: developing a visual diagnostic interface that synchronously renders and displays dynamic stiffness spectrum, modal vibration animation and stress cloud map of key areas. Users can select a specific frequency point through an interactive slider to view the magnification of local deformation, which can be magnified up to 50 times to identify weak links. The interface also provides a stiffness attenuation rate index, defined as the percentage decrease of the average dynamic stiffness in a certain frequency band relative to the benchmark value. When it exceeds 15%, a warning indicator is triggered.

[0014] Preferably, the method is applied in the lightweight design iteration process of the reducer housing. After each topology optimization or size optimization, this calculation process is automatically called to perform dynamic performance verification, forming a closed-loop optimization system. The optimization objective function includes two indicators: mass reduction rate and dynamic stiffness retention rate under the third bending mode. The constraints cover the maximum static stress, fatigue safety factor and the lowest first natural frequency, ensuring that the structure meets the dynamic stiffness requirements while reducing weight.

[0015] Preferably, it also includes: setting up a multi-level verification mechanism. The first level uses laboratory bench sweep frequency test to obtain measured dynamic stiffness data, arranges no less than 8 triaxial acceleration sensors in the area around the bearing housing, applies a sinusoidal signal with a sweep frequency range of 10 Hz to 2 kHz to the exciter, collects the transfer function and compares it with the simulation results, with a consistency correlation coefficient of no less than 0.88; the second level introduces a digital twin system to collect the vibration signal of the reducer housing in real time during the vehicle durability test, and performs dynamic stiffness tracking by executing a simplified version of the method online through edge computing devices to realize performance degradation monitoring under service conditions.

[0016] Compared with the prior art, the present invention has the following beneficial effects:

[0017] To improve the accuracy of dynamic stiffness calculation, the introduction of manufacturing process-related residual stress fields and material anisotropy characteristics into the finite element model significantly improves the stiffness overestimation problem caused by traditional equivalent boundary conditions, increasing the correlation between simulation results and measured data from less than 0.7 to over 0.9. The multiphysics coupled excitation model fully covers the real load environment in the operation of electric drive systems, solving the dynamic response distortion defects caused by single load assumptions, especially reducing the phase lag prediction error by more than 40% in the high-frequency band above 1 kHz.

[0018] The integrated dynamic analysis framework abandons the traditional process of separating modal analysis and harmonic response analysis, and adopts a direct solution method of nonlinear transient dynamics. It fully preserves the nonlinear characteristics and time evolution information of the system, and realizes end-to-end calculation from the original excitation to the extraction of complex stiffness parameters, avoiding information loss in intermediate links. The extraction of complex stiffness parameters not only gives the amplitude response, but also accurately characterizes the energy dissipation mechanism, providing a theoretical basis for the layout of damping materials and structural modification.

[0019] Enhanced engineering application adaptability, a quantitative mapping relationship between process parameters and structural dynamic performance was established, enabling the impact of manufacturing deviations on dynamic stiffness to be predicted during the design phase, supporting the implementation of manufacturing-oriented design concepts; the integration of closed-loop optimization system and visual diagnostic tools significantly shortens the development cycle, with a single analysis iteration time controlled within 40 minutes, improving efficiency by more than 6 times compared to the traditional trial-and-error mode, and is suitable for the rapid R&D needs of new energy vehicle reducers.

[0020] It supports full lifecycle performance management. By embedding the online computing capabilities of the digital twin system, it enables dynamic tracking of the dynamic stiffness of the reducer housing under actual working conditions, breaking through the limitation that traditional offline analysis cannot reflect performance degradation. The early warning mechanism can identify stiffness deterioration trends caused by fatigue damage or loose connections in advance, providing decision support for preventive maintenance and extending product service life. Attached Figure Description

[0021] Figure 1 This is a schematic diagram of the overall technical solution architecture of a method for calculating the dynamic stiffness of the bearing seat of a new energy reducer proposed in this invention. Detailed Implementation

[0022] Please refer to Figure 1 To make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to specific embodiments.

[0023] In the above method for calculating the dynamic stiffness of the bearing seat of the new energy reducer housing, step (1) constructs a refined finite element model. Its goal is to solve the problem of structural response prediction distortion caused by neglecting the directionality of material microstructure and manufacturing process disturbances in traditional modeling methods, thereby providing a high-fidelity numerical carrier for subsequent dynamic performance analysis. Specifically, step (1) is based on the three-dimensional geometric model of the new energy reducer housing. This three-dimensional geometric model comes from the CAD design file, in STEP or IGES format, with an accuracy level not lower than ISO 10303-21 standard, and contains complete topological information and geometric dimension parameters. Based on this model, a high-density tetrahedral and hexahedral hybrid mesh generation strategy is implemented. Hexahedral elements are preferentially used for regular geometric regions such as the bearing seat body and the main wall of the housing to ensure calculation accuracy and convergence; tetrahedral elements are used to handle complex curved surface regions such as stiffener intersection areas and casting fillet transition zones to ensure geometric conformity. The overall mesh density is controlled at no less than 150 nodes per cubic millimeter, and the minimum element size is strictly controlled within 0.5 mm. Local densification is implemented, especially in key locations such as the inner ring contact area of ​​the bearing housing and the bolt preload application point, so that the element side length is reduced to the range of 0.3 mm to 0.4 mm, in order to fully capture stress gradient changes.

[0024] An adaptive meshing strategy based on geometric curvature and gradient rate of change was adopted during the mesh generation process. This strategy calculates the normal curvature value of each micro-element on the model surface in real time. When a curvature greater than 0.02 mm is detected in a certain region, a local refinement mechanism is automatically triggered, reducing the mesh size of that region to 60% of the original base mesh, and transitioning outward layer by layer to ensure that the size ratio of adjacent elements does not exceed 1.2 times, avoiding numerical ill-conditioned phenomena caused by abrupt size changes. Simultaneously, multiple mesh refinement layers are set at the edges of bolt connection holes, the junction of stiffeners and the box body, and areas with casting fillet radii less than 3 mm. Each layer has a thickness decreasing by 15%, forming a progressive transition structure that effectively alleviates the discretization error caused by stress concentration. Mesh quality evaluation indicators include a Jacobi ratio of not less than 0.7, a warpage of less than 15°, and an aspect ratio controlled within 1:5. All unqualified elements are corrected using a local re-meshing algorithm.

[0025] In the material property definition stage, an anisotropic elastic modulus tensor is introduced to characterize the microstructure orientation of die-cast aluminum alloys. This tensor takes the form:

[0026]

[0027] in, For the anisotropic elastic modulus tensor, It is the axial tensile stiffness along the x-direction (i.e., the ratio of normal stress to strain generated when stretched in the x-direction). Axial tensile stiffness along the y-direction , and This is a term related to Poisson's ratio, representing the degree of lateral contraction in one direction when the object is stretched in another (e.g., ...). (This represents the stiffness that is compressed in the y direction due to tension in the x direction). This represents the axial tensile stiffness along the z-direction (i.e., perpendicular to the die-casting flow direction), reflecting the material's tensile strength in the vertical direction. This refers to the shear stiffness (shear modulus in the yz plane). This refers to shear stiffness (shear modulus in the xz plane). The shear stiffness (shear modulus in the xy plane) is calibrated based on measured data of cooling rate and mold orientation during the die-casting process. The value range is: 78.5 GPa to 82.3 GPa. The value range is: 26.4 GPa to 28.7 GPa. The value range is: 29.1 GPa to 31.5 GPa. The value range is: 26.8 GPa to 28.9 GPa. The value range is: 27.2 GPa to 29.4 GPa. The value ranges from 26.1 GPa to 28.3 GPa, reflecting the difference in mechanical properties along the casting direction and the perpendicular direction. Furthermore, to describe the plastic yielding behavior of the material under cyclic loading, the Chaboche viscoplastic constitutive model is introduced, whose expression is:

[0028]

[0029] in, For plastic strain rate, The initial yield strength, For equivalent plastic strain rate, To accumulate equivalent plastic strain, It is a nonlinear hardening modulus. The back stress evolution rate, For the deviatoric stress tensor, For the back stress tensor, a total of 3 independent sets were set ( , The parameter combinations correspond to the hardening behavior under three working conditions: low-cycle fatigue, mid-frequency vibration, and high-frequency impact, respectively, and the nonlinear hardening modulus. The value range is from 850 MPa to 1100 MPa. The values ​​range from 80 to 120, determined through inversion from uniaxial tensile-compression cyclic test data. This model accurately simulates the Bauschinger effect and ratchet strain accumulation process, ensuring good predictive consistency even under non-proportional loading paths. The final generated finite element model has a total of 1.2 million to 1.8 million nodes, with element types covering C3D8R (eight-node linear hexahedron), C3D10M (ten-node quadratic tetrahedron), and S4R (four-node shell element). Pre-processing was performed using the ABAQUS CAE platform, outputting INP format input files for the solver to use.

[0030] In the above method for calculating the dynamic stiffness of the bearing seat of the new energy reducer housing, step (2) integrates the manufacturing process disturbance field. Its goal is to transform the uncontrollable process variables in the actual production process into quantifiable physical field inputs, thereby improving the realism of the simulation boundary conditions. Specifically, step (2) first obtains the historical process parameter data stream accumulated on the die-casting production line of the new energy reducer housing. This data stream comes from the MES system database, with a sampling frequency of 1 time / second, and a time span covering no less than 500 batches of production records in the most recent 12 months. Each batch includes four core variables: mold temperature, pouring speed, holding time, and cooling rate. The mold temperature monitoring points are arranged at 5 key locations on the cavity surface, using K-type thermocouple sensors with a measurement range of 150℃ to 300℃ and an accuracy of ±1.5℃. The pouring speed is obtained from feedback by the servo control system, in meters per second, with a resolution of 0.01 m / s. The holding time is recorded with an accuracy of 0.1 seconds. The cooling rate is obtained by capturing the temperature decay curve within 10 seconds before and after demolding using an infrared thermometer array, and then calculating it by differentiation, in degrees Celsius per second.

[0031] Based on the aforementioned historical data, a residual stress prediction neural network model was constructed. This model employs a ConvLSTM hybrid architecture with a depth of 6. The first three layers are one-dimensional convolutional layers with a kernel size of 5 and a stride of 1, using ReLU activation to extract spatial correlation features of different process parameter sequences. The last three layers are LSTM layers with 128 hidden units. The gating mechanism includes forget gates, input gates, and output gates to capture the dependencies and dynamic evolution patterns in the time series. The training set consists of no fewer than 50,000 samples, each containing a 60-second time window (i.e., the process parameter sequence for the first 60 seconds) and corresponding residual stress distribution labels. The label data comes from the actual measurement results of the same batch of shells on an X-ray diffraction device. The measurement locations cover nine typical cross-sections, including the outer ring of the bearing housing, the flange connection area, and the bottom support foot. Each cross-section has no fewer than six measurement points, with a measurement depth of 0.1 mm. The instrument model is Bruker D8 ADVANCE, and the measurement error is controlled within ±10 MPa. The training process uses the Adam optimizer with an initial learning rate of 0.001, a batch size of 64, and a loss function of mean squared error (MSE). After at least 200 iterations, the model's coefficient of determination on the validation set is calculated. It remains stable above 0.93, meeting the requirements for engineering applications.

[0032] After model training, it is deployed on a high-performance computing server to receive real-time input of the current batch of process parameters and output a three-dimensional residual stress distribution field inside the shell. This output is stored in voxelized data format with a spatial resolution of 2 mm × 2 mm × 2 mm, containing approximately 150,000 voxel elements. Each element records three normal stress components (σxx, σyy, σzz) and three shear stress components (τxy, τyz, τzx). To achieve field quantity mapping with the finite element model, a spatial coordinate alignment operation is performed, aligning the origin of the voxel mesh with the origin of the global coordinate system of the finite element model. A trilinear interpolation algorithm is then used to map the stress values ​​at the voxel nodes to the nodes of the finite element solid mesh, with interpolation weights calculated based on the inverse distance weighting principle. To eliminate field quantity jumps caused by discretization, moving least squares (MLS) is further used for smoothing. The search radius is set to 5 mm, and a quadratic polynomial is used as the basis function. After smoothing, the standard deviation of the residual stress field is reduced by more than 35%. Finally, the predicted residual stress field is superimposed as a prestress field into the finite element model established in step (1), and imported through the "Initial Conditions" module to form an initial state containing manufacturing disturbances. The calibration process adopts an iterative feedback mechanism: three representative working conditions are selected for simulation and actual measurement comparison. If the maximum residual stress deviation exceeds 8%, the process returns to adjust the hyperparameters of the neural network and retrains until the error control requirements are met.

[0033] In the above method for calculating the dynamic stiffness of the bearing seat of the new energy reducer housing, step (3) couples multi-physics dynamic loads. Its goal is to reproduce the comprehensive excitation of the new energy electric drive system under real operating conditions and avoid the distortion of dynamic response caused by single load assumptions. Specifically, step (3) constructs a comprehensive excitation model under the operating conditions of the electric drive system, and simultaneously applies three main types of external excitations: the first type is the periodic alternating torque caused by gear meshing, which originates from the time-varying meshing stiffness and transmission error of the gear pair between the high-speed and low-speed stages of the reducer, and is obtained through simulation using Romax Designer software. The excitation frequency is the main meshing frequency and its harmonics (f_m = z_1 × n_1 / 60, where z_1 is the number of teeth of the pinion and n_1 is the rotational speed), with an amplitude range of 120 N·m to 480 N·m, and the phase difference is set according to the gear modification parameters; the second type is the high-frequency axial ripple load generated by the electromagnetic pulsation force of the motor, which is modeled based on the harmonic components of the stator winding current of the motor, mainly considering the 6th, 12th, and 18th spatial harmonic components, whose amplitudes are obtained by inversion using the electromagnetic simulation software JMAG, and are 85 N, 42 N, and 28 N respectively under rated operating conditions. N, the phase difference is dynamically adjusted based on feedback information from the rotor position sensor, with an update cycle of 0.5 milliseconds to ensure synchronization with the rotation angle; the third type is the random vibration displacement boundary condition where road excitation is transmitted to the suspension point during vehicle movement. This excitation is generated based on the Class C road power spectral density function (PSD), and its expression is:

[0034]

[0035] Among them, the reference space frequency =0.1 Reference value

[0036] =256× cubic meters, index =2, the vehicle speed is set to 60 km / h, a time-domain random signal is generated by inverse Fourier transform, and then the transmission path is reversed through the multibody dynamics model of the whole vehicle established by Adams / Car to obtain the six-degree-of-freedom displacement input of the four suspension points of the reducer. The sampling frequency is 5 kHz and the duration is 60 seconds.

[0037] The above three types of loads are synchronously applied in the Abaqus / Explicit environment. The gear meshing torque is applied to the input shaft end face as a torque load, with the loading direction tangential to the axis. The electromagnetic pulsating force is converted into an equivalent harmonic force acting on the central axis of the bearing housing, with the direction axial. Dynamic amplitude and phase updates are implemented through DLOAD subroutine programming. The road surface excitation is applied to the suspension mounting hole node under BOUNDARY conditions, defined as a displacement function varying with time, and imported using a .tab file. The temperature field is imported through the thermo-mechanical coupling analysis module, with the working temperature rise range set to 25℃ to 120℃. The heating process follows a linear law, taking 300 seconds. The material's thermal expansion coefficient varies linearly with temperature in a piecewise manner, reaching 21.5 × 10⁻⁶ in the 25℃ to 60℃ range. / ℃, in the range of 60℃ to 100℃, it is 22.8× / ℃, in the range of 100℃ to 120℃, it is 23.6×10 / ℃, all physical properties are defined via MATERIAL DATA TABLE. Boundary constraints are set as follows: release the corresponding degrees of freedom at the suspension pivot to match the actual rubber bushing characteristics, and apply full-degree-of-freedom fixed constraints to the remaining mounting holes. The total loading process time is set to 1.2 seconds, covering at least 5 complete road excitation cycles and hundreds of gear meshing events to ensure the system enters a steady-state response.

[0038] In the above method for calculating the dynamic stiffness of the bearing seat of the new energy reducer housing, step (4) performs a nonlinear transient dynamic solution. Its goal is to capture the full-time response of the system under strong nonlinearity and multi-excitation coupling, and retain complete dynamic evolution information. Specifically, step (4) adopts the explicit central difference integral algorithm as the basic time-progression framework, and draws on the high-frequency damping control idea of ​​the Hilber-Hughes-Taylor (HHT) method to introduce a numerical dissipation mechanism to suppress non-physical high-frequency oscillations caused by mesh discretization or nonlinear contact. The explicit integral algorithm uses an improved version of the central difference method, and its basic recursive formula is:

[0039]

[0040]

[0041] in, , , The first The displacement, velocity, and acceleration vectors of each step. For the time step, it should be noted that classic HHT methods are typically based on implicit Newmark formats, using parameters... , , Jointly define the integral weights, where , However, in large-scale explicit simulation scenarios, implicit iteration significantly increases computational overhead. Therefore, this invention does not directly adopt the standard HHT implicit scheme, but instead retains the efficient explicit structure, only referencing it through a single parameter. The core mechanism for controlling high-frequency dissipation.

[0042] In terms of specific implementation, in completing acceleration After explicit calculation, an HHT-type damping correction is introduced:

[0043]

[0044] in, The effective acceleration after damping correction. The numerical damping coefficient specified by the user, with a value range of [value range missing]. 0.05≤ ≤0. When When = 0, it degenerates into the undamped central difference method; when When <0, high-frequency modes are selectively attenuated, while low-frequency dynamic behavior is basically unaffected.

[0045] Time step An adaptive control strategy is adopted, with an initial value of 1 μs, and strictly satisfies the Courant-Friedrichs-Lewy (CFL) stability condition:

[0046]

[0047] in, The characteristic length (in meters) of the smallest element in the finite element model. The longitudinal wave velocity in the shell material (typically die-cast aluminum alloy) is approximately 5100 m / s. In actual calculations, the minimum time step can be dynamically reduced to 0.2 μs, corresponding to an equivalent sampling frequency of no less than 5 MHz, thereby ensuring accurate resolution of key dynamic response components up to 10 kHz.

[0048] In terms of solver configuration, double-precision floating-point operation mode is enabled, with a memory allocation of no less than 128 GB. The MPI parallel computing framework is used, divided into 8 computational domains. The communication protocol is InfiniBand with a bandwidth of 40 Gbps, and the load balancing efficiency is higher than 92%. An energy conservation monitoring module is simultaneously set up during the solution process. This module calculates the total system energy E_total = E_kinetic + E_internal + E_external every 100 steps, where the kinetic energy E_kinetic is calculated from the sum of squared nodal velocities, the internal energy E_internal is obtained from the element stress-strain integral, and the external work E_external is accumulated from the work done by the loads. When the fluctuation of E_total relative to the initial energy exceeds 3%, it is considered a risk of energy divergence. The system automatically reduces the current time step to 70% of the original value and restarts the local calculation within that time window, allowing a maximum of 3 consecutive restarts. If convergence is still not achieved, the condition is marked as abnormal and the solution is terminated. All computational tasks are submitted to the Linux cluster scheduling system Slurm, with the average single solution time controlled within 38 minutes. The final output includes: displacement response sequences and corresponding reaction force time history curves of no less than 200 key nodes on the bearing housing mounting surface. The data format is CSV, with no less than 240,000 sampling points. It contains complete information on the three degrees of freedom: radial (X direction), axial (Y direction), and tangential (Z direction). The timestamp accuracy reaches the microsecond level, providing the original data foundation for subsequent frequency domain analysis.

[0049] In the above method for calculating the dynamic stiffness of the bearing seat of the new energy reducer housing, step (5) extracts complex dynamic stiffness parameters. Its goal is to separate dynamic stiffness characteristics with clear physical meaning from the complex time-domain response, supporting structural performance evaluation and optimization decisions. Specifically, step (5) first performs synchronous Fourier transform on the displacement and reaction force signals obtained in step (4) to obtain the transfer function matrix H(ω) in the frequency domain, where the element H_ij(ω) represents the response amplitude and phase of the i-th degree of freedom when a unit excitation is applied to the j-th degree of freedom. To prevent spectral leakage, the original signal is weighted by a Hanning window before the transformation. The window function length is equal to the total signal duration, the overlap rate is 50%, and the frequency resolution is controlled within 1 Hz to ensure that adjacent modal peaks can be clearly distinguished. The transformed data is stored as a complex matrix with a dimension of 3N×3N (N is the number of key nodes), and batch processing is completed using the MATLAB Signal Processing Toolbox.

[0050] Subsequently, modal parameters were fitted to the transfer function matrix using the least squares complex frequency domain method (LSCF), limiting the modal order to the range of 1 to 15, covering the key frequency band from 0 Hz to 2000 Hz. A weighted least squares criterion was used during the fitting process, with the weight function W(f) assigned according to the signal-to-noise ratio (SNR) distribution. The weight was set to 0.3 for frequency bands with an SNR below 20 dB, 0.6 for 20–40 dB, and 1.0 for frequencies above 40 dB, prioritizing the fitting accuracy in the formant region. The fitting output includes the natural frequencies of each mode. (Unit: Hz), Damping Ratio (Dimensionless), mode participation factor (Unit: mm / N) and complex mode vector Based on this, the complex dynamic stiffness is further calculated. The real part Represents the energy storage stiffness, imaginary part The energy consumption characteristics are characterized by the following formula:

[0051]

[0052] in, Angular frequency, These are the modal orders, from order 1 to order 15. For the first Mode participation factor of first-order modes For the first The transpose of the first-order mode shape vector. Represents transpose. For the first Damping ratio of first mode, For the first The natural angular frequency of the first mode, i.e. The complex stiffness matrix was normalized according to the international standard ISO 10846, converting the original unit N / m to kilonewtons per micrometer (kN / μm) for easier comparison across projects. Two types of visualizations were generated: one is a frequency-dynamic stiffness spectrum, with the horizontal axis representing frequency (10 Hz to 2000 Hz) and the vertical axis representing the complex stiffness amplitude |K*|, displayed on logarithmic coordinates, marking the resonant frequency points of each order; the other is a Nyquist plot, plotting the complex plane trajectory with K″ as the vertical axis and K′ as the horizontal axis, where each closed-loop curve corresponds to a first-order mode, used to determine the system's stability and damping characteristics. All charts were generated using the Python Matplotlib library, with a resolution of at least 300 dpi, and saved in both PDF and PNG formats.

[0053] The method described above for calculating the dynamic stiffness of the bearing seat in the new energy reducer housing also includes establishing a process-structure-performance correlation database, aiming to solidify knowledge assets and support rapid prediction. This database is built using MySQL 8.0 with a utf8mb4 character set and an InnoDB storage engine. It contains the following core table structures: the "materials_params" table records anisotropic elastic modulus, Chaboche model parameters, etc., with fields including C11, C12…γ3, and a data type of DOUBLE; the "process_conditions" table records process variables such as mold temperature and casting speed, with fields such as temperature and speed, and a data type of FLOAT; the "mesh_config" table records the total number of meshes, minimum element size, and coordinates of the encrypted area, with data types of INT and JSON; the "dynamic_stiffness_curve" table stores complex dynamic stiffness frequency response data, using a BLOB type to store serialized NumPy arrays, supplemented by JSON fields indicating the frequency range and sampling interval. The database is incrementally backed up daily to an off-site storage server, with a retention period of no less than 5 years. The structured dataset generated from this database was used to train a regression-type support vector machine (SVM) model. The kernel function used was radial basis function (RBF), with a penalty coefficient C=100 and γ=0.01. The training set size was no less than 8000 sets, and the test set accounted for 20%. After training, the model was deployed on a Web API service to receive input parameters from new design schemes. The prediction response time was less than 3 seconds, and the relative error was less than 7%, significantly accelerating the performance evaluation in the early conceptual design stage.

[0054] The method for calculating the dynamic stiffness of the bearing housing of the new energy reducer also includes the development of a visual diagnostic interface, aiming to improve human-computer interaction efficiency and defect identification capabilities. This interface is developed based on the Qt 5.15 framework, written in C++, and uses OpenGL 4.5 as its graphics rendering engine, supporting both Windows and Linux platforms. The main view area of ​​the interface simultaneously renders three items: the left side displays a frequency-dynamic stiffness spectrum, allowing users to select a specific frequency point (accuracy 0.1 Hz) using an interactive slider; the upper right side displays a modal shape animation, showing the overall deformation of the housing at that frequency, with the deformation magnification continuously adjustable from 1x to 50x to identify local weak points; and the lower right side displays a stress cloud map of the key area, with the color mapping range automatically adjusted according to the maximum von Mises stress, 11 color levels, and support for contour overlay display. The bottom panel of the interface provides a function to calculate the stiffness attenuation rate index, defined as the percentage decrease in the average dynamic stiffness within a certain frequency band (e.g., 800 Hz to 1200 Hz) relative to a reference value (e.g., stiffness at 100 Hz). The calculation formula is as follows:

[0055]

[0056] when When the error exceeds 15%, the system automatically adds a flashing red warning indicator to the spectrum and pops up a prompt box listing possible causes (such as insufficient local wall thickness, unreasonable stiffener layout, etc.). The interface supports parallel loading and comparative analysis of multiple projects, and all operation logs are recorded to the system log file for easy auditing and tracing.

[0057] In the above method for calculating the dynamic stiffness of the bearing seat in a new energy reducer housing, the method is applied to the iterative process of lightweight design for the reducer housing, aiming to achieve optimal weight reduction under dynamic performance constraints. Specifically, after each topology optimization or dimensional optimization, this calculation process is automatically invoked for dynamic performance verification, forming a closed-loop optimization system. The optimization process is driven by the Isight platform, with upstream connection to OptiStruct for structural optimization and downstream connection to Abaqus to execute the method described in this invention. The optimization objective function is defined as:

[0058]

[0059] in, For the current quality, For the initial mass, This represents the dynamic stiffness retention rate under the third bending mode. As the benchmark value, the weighting coefficient =0.6, =0.4. Constraints include: maximum static stress not exceeding 70% of the material's yield strength (i.e., ≤180 MPa), fatigue safety factor not less than 1.8 (calculated based on Miner's linear cumulative damage theory), minimum first-order natural frequency not less than 850 Hz, and prevention of resonance with the excitation source. After each optimization iteration, the system automatically extracts a new geometric model, re-meshes the grid, and performs a full-process calculation to determine if all constraints are met. If not, it returns to adjust design variables until a Pareto optimal solution is found. The entire closed-loop process averages 12 iterations, with each analysis iteration taking less than 40 minutes, improving efficiency by more than 6 times compared to traditional trial-and-error methods.

[0060] The method for calculating the dynamic stiffness of the bearing housing of a new energy reducer also includes a multi-level verification mechanism to ensure the reliability and engineering applicability of the simulation results. The first level of verification employs a laboratory bench sweep frequency test conducted on an electric vibration table. The exciter is an LDS V994 with a maximum thrust of 10 kN. The sweep frequency range is 10 Hz to 2000 Hz, and the acceleration amplitude is controlled within 5 g to avoid nonlinear interference. The housing under test is fixed to the table via four suspension points. At least eight triaxial accelerometers (PCB 356A16, sensitivity 100 mV / g) are evenly distributed around the bearing housing, with a sampling frequency of 5 kHz. The transfer function is acquired using LMS Test.Lab software. Simulation results After undergoing the same filtering process, the data is compared with the measured data to calculate the consistency correlation coefficient. :

[0061]

[0062] in, As a modal guarantee criterion, To simulate the frequency response function, This is the measured frequency response function. This is the Hermitian transpose (conjugate transpose). To simulate the autocorrelation term, To measure the autocorrelation term, the following is required: The value must be at least 0.88; otherwise, the model parameters must be corrected. The second level of verification introduces a digital twin system. During the vehicle durability test, a wireless vibration sensor (sampling rate 2 kHz) installed on the reducer housing collects vibration signals in real time. The data is transmitted via 4G / 5G network to an edge computing device (NVIDIA Jetson AGX Xavier), where a simplified method (reduced-order model ROM + fast FFT) is executed online. The dynamic stiffness tracking curve is updated every 5 minutes to monitor performance degradation under service conditions. When a stiffness drop in a certain frequency band exceeds a threshold (e.g., 15%) or a new resonance peak appears, a cloud-based alert is triggered, notifying maintenance personnel to arrange inspections, thereby extending product lifespan.

[0063] Building upon this foundation, a closed-loop health management platform based on an edge-cloud collaborative architecture is further constructed to achieve fully automated response across the entire chain, from data acquisition and model inference to decision support. This platform consists of three parts: a front-end perception layer deployed on a high-precision MEMS vibration sensor array on the test vehicle, with a sampling frequency of 4 kHz and self-calibration capabilities, capable of real-time compensation for temperature drift and installation angle deviations; a middle computation layer employing a lightweight reduced-order model (ROM), compressing the original million-degree-of-freedom finite element model to a state-space expression of no more than order 50 through modal truncation and Krylov subspace projection methods, ensuring that a single dynamic stiffness calculation on the Jetson AGX Xavier edge device takes less than 800 milliseconds, meeting near real-time processing requirements; and a back-end analysis layer located on a private cloud server, receiving aggregated data from multiple test vehicles and performing long-term trend analysis and group failure mode mining.

[0064] During online operation, the edge device performs a sliding window FFT analysis every 30 seconds to extract the main vibration energy distribution frequency band of the bearing housing area under the current operating conditions. It then performs dynamic comparison with the pre-stored health benchmark spectrum to calculate the stiffness attenuation rate index η(t). If three consecutive measurement results show that η(t) > 15%, a deep diagnostic process is initiated: First, the ROM model is called to inject the current measured excitation conditions (including CAN bus signals such as motor speed, torque, and vehicle speed) to simulate the structural response and invert the equivalent stiffness parameters. Then, using a transfer learning strategy, the residual stress-stiffness sensitivity matrix calibrated in the laboratory is mapped to the current service environment to preliminarily determine the possible failure mechanism. For example, if the stiffness in the low-frequency band (<500 Hz) decreases significantly and is accompanied by local temperature rise, it is likely to be judged as bolt preload loosening. If the resonance peak in the high-frequency band (>1.5 kHz) shifts significantly, it suggests that there may be microcrack propagation caused by casting defects.

[0065] All diagnostic results are accompanied by a confidence score (derived from historical false alarm rates) and uploaded to a cloud database via an encrypted communication protocol. The cloud platform utilizes a big data analytics engine to perform horizontal comparisons of data from multiple vehicles and mileage segments to identify regionally common problems. For example, in a batch of vehicles undergoing high-altitude durability testing, samples at altitudes above 3000 meters generally exhibited earlier stiffness degradation, which, after correlating with meteorological data, was inferred to be related to accelerated material fatigue due to large diurnal temperature variations. Such insights are fed back to the design department to optimize thermal-mechanical coupling protection measures for subsequent products.

[0066] Furthermore, the system integrates a knowledge graph module, establishing semantic associations between fault cases, maintenance records, simulation parameters, and sensor data, supporting natural language queries. Engineers can input "find cases of abnormal noise caused by decreased bearing housing stiffness within the last three months" via voice input, and the system automatically matches similar vibration characteristics with maintenance reports to assist in quickly locating the root cause. The entire digital twin system not only enhances product lifecycle management capabilities but also provides a valuable empirical feedback loop for the calculation method described in this invention, driving continuous iterative optimization of the simulation model.

[0067] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for calculating the dynamic stiffness of a new energy retarder housing bearing seat, characterized in that, The method comprises the following steps: A refined finite element model is constructed based on a three-dimensional geometric model of the new energy decelerator shell, and an anisotropic elastic modulus tensor and a nonlinear constitutive relationship are introduced in the material properties; A manufacturing process disturbance field is integrated, which is a three-dimensional residual stress distribution field predicted by a neural network model trained based on pressure casting process parameter historical data flow, and is superimposed on the refined finite element model as a prestress field; Multi-physical field dynamic loads are coupled, including periodic alternating torque caused by gear meshing, high-frequency axial fluctuating load generated by motor electromagnetic pulsation force, and random vibration displacement boundary conditions of the suspension point conducted from the road excitation in the vehicle driving, and a temperature field is synchronously introduced; Nonlinear transient dynamics solving is performed to obtain displacement response sequences of each node on the bearing seat mounting surface in multiple degrees of freedom and corresponding time history curves of the counterforce; The displacement response sequences and the counterforce time history curves are processed to extract complex dynamic stiffness parameters.

2. The method of claim 1, wherein, The refined finite element model is constructed, comprising: The three-dimensional geometric model is divided into high-density tetrahedral and hexahedral mixed grids, and local grid densification is implemented in the bearing seat area, bolt connection hole, rib plate transition zone and casting fillet position; In the material property definition, the anisotropic elastic modulus tensor is used to represent the microstructure directionality of the die-cast aluminum alloy, and the viscoplastic constitutive model containing multiple sets of back stress tensors is used to describe the cyclic plastic behavior of the material.

3. The method of claim 2, wherein, The high-density tetrahedral and hexahedral mixed grids are divided, comprising: An adaptive division strategy based on geometric curvature and gradient change rate is adopted to automatically trigger local densification in areas with curvature exceeding a threshold value, and the size ratio of adjacent cells is controlled to avoid numerical ill-conditioning.

4. The method of claim 1, wherein, The manufacturing process disturbance field is integrated, comprising: Mold temperature, pouring speed, holding time and cooling rate are obtained as input variables and input into the neural network model; The voxelized residual stress field output by the neural network model is mapped to the nodes of the refined finite element model through spatial coordinate alignment and interpolation algorithm, and is smoothed to eliminate field quantity jumps.

5. The method of claim 1, wherein, The multi-physical field dynamic loads are coupled, comprising: The motor electromagnetic pulsation force is modeled according to the harmonic components of the motor stator winding current, and its phase is dynamically adjusted according to the rotor position information; A random signal is generated based on the road power spectral density function, and a transfer path is back calculated via a vehicle multibody dynamics model to form the displacement boundary conditions acting on the decelerator suspension support point.

6. The method of claim 1, wherein, The nonlinear transient dynamics solving is performed, comprising: An explicit integration algorithm combined with damping correction technique is used for full-time domain response calculation; An energy conservation monitoring module is set, and when the total energy fluctuation of the system exceeds a preset threshold, the time step is automatically reduced and the local calculation is restarted.

7. The method of claim 1, wherein, The displacement response sequences and the counterforce time history curves are processed to extract complex dynamic stiffness parameters, comprising: Synchronous Fourier transform is performed on the displacement response sequences and the counterforce time history curves to obtain a transfer function matrix in the frequency domain; Modal parameter fitting is performed on the transfer function matrix based on the least squares complex frequency domain method to obtain natural frequency, damping ratio and mode shape participation factor; According to the modal parameters, complex stiffness is calculated, and a real part of the complex stiffness represents a storage stiffness, and an imaginary part represents a dissipation characteristic.

8. The method of claim 7, wherein, Synchronous Fourier transform is performed on the displacement response sequence and the time history curve of the counterforce, including: Before the transform, a window function weighting process is performed on the original signal to prevent spectrum leakage and control frequency resolution.

9. The method of claim 1, wherein, Further comprising: A process-structure-property correlation database is established to store material parameters, process conditions, mesh configurations, and dynamic stiffness curves; Based on the process-structure-property correlation database, a machine learning prediction model is trained to quickly predict the dynamic stiffness frequency response characteristics of a new design scheme.

10. The method of claim 1, wherein, Further comprising: A visual diagnosis interface is developed to synchronously render dynamic stiffness spectrograms, modal vibration mode animations, and stress cloud maps of key areas; In the visual diagnosis interface, an interactive tool is provided for a user to select a specific frequency point to view local deformation, and a stiffness attenuation rate index is calculated to trigger an early warning.

Images