A monte carlo prediction method for assembly accuracy considering micro-topography and contact deformation

Abstract

The application relates to the field of precision mechanical assembly and tolerance analysis, and discloses a kind of assembly precision Monte Carlo prediction method considering micro topography and contact deformation. For inertial platform ring-column assembly interface, small displacement rotation macroscopic pose deviation parameters are established, and micro topography deviation parameters based on Zernike and Legendre-Fourier polynomials are established. The contact deformation under pre-tightening force is solved by conjugate gradient-fast Fourier transform method, and the axis position, inclination and coaxiality index are fitted. A thousand times of Monte Carlo simulation is performed by using a multi-process parallel framework to generate a sample data set. The parameter contribution proportion is quantified by variance contribution decomposition, and the sensitivity ranking is output. The application breaks through the limitation of traditional rigid body analysis, reduces the calculation time to seconds, accurately identifies the key tolerance elements, and provides quantitative basis for form and position tolerance design and process optimization.

Description

Technical Field

[0001] This invention relates to the field of precision mechanical assembly and tolerance analysis technology, specifically a Monte Carlo method for predicting assembly accuracy that takes into account microstructure and contact deformation. Background Technology

[0002] Precision inertial platform systems are core components in aerospace navigation and guidance equipment. A typical structure consists of a central platform, inner and outer sun rings, and multiple pairs of intermediate axes, all connected and fixed by bolts. During system operation, the central platform and the inner and outer sun rings rotate at high speed relative to each other along their respective axes. The orthogonality accuracy between these two rotation axes directly affects the system's attitude calculation accuracy and navigation performance. Since axes are not measurable physical features, their spatial positions are entirely determined by the actual contact state of the assembly interfaces. Therefore, assembly accuracy is a key technical factor determining the product's functional specifications in such systems. Assembly accuracy analysis involves the coupling and transmission of multi-source deviations in engineering. Part dimensions fluctuate randomly within tolerance limits, mating surfaces exhibit microscopic machining textures, bolt tightening causes elastic deformation of the contact surfaces, and errors between multiple assembly interfaces can be transmitted and accumulated across interfaces. How to predict the assembly accuracy under the combined effects of these factors during the design phase is a core research problem in the field of precision mechanical assembly.

[0003] Currently, there are three main types of methods for analyzing assembly accuracy. The first type is the extreme value method, which assumes that all deviations are simultaneously calculated by superimposing their most unfavorable extreme values. While simple, this method yields overly conservative results, leading to tolerances far exceeding actual requirements and causing unnecessary increases in processing costs. The second type is the statistical tolerance method, represented by the root mean square method. This method assumes that each error source is independent and follows a normal distribution, estimating assembly accuracy through root mean square synthesis. It is computationally efficient but struggles to handle nonlinear error propagation and non-normally distributed inputs. The third type is the deterministic analysis method based on finite element simulation. This method can simulate contact deformation and stress distribution under specific assembly conditions, but each simulation is time-consuming and cannot statistically cover a large number of possible state combinations within the tolerance range.

[0004] The existing patent, CN121960017A, is titled "A Method for Analyzing Assembly Accuracy of Rotary Shaft System of Mechanical Structure Based on Laser Tracking Measurement System." This method uses the Jacobi-spinometer model to analyze assembly accuracy. It establishes a local coordinate system on each mating surface of the assembly, uses the SDT six-parameter model to describe the spatial pose deviation of each mating surface, and linearly superimposes each deviation along the transfer chain to the functional target using the Jacobi matrix. Finally, it uses the Monte Carlo method to statistically analyze the distribution range of assembly accuracy.

[0005] Existing assembly accuracy analysis methods have the following four shortcomings.

[0006] First, existing methods do not fully characterize the microstructure of mating surfaces. Machined surfaces are not perfectly smooth at the microscale; instead, they exhibit a wave-and-trough structure caused by the machining process. Under bolt preload, the micro-protrusions on the two mating surfaces undergo elastic contact deformation, and the number, area, and distribution of contact spots directly affect the positioning accuracy of the axis after assembly. Existing methods either ignore surface morphology or describe it using only a single roughness parameter, failing to establish a complete computational model from surface microstructure to contact state and then to axis positioning.

[0007] Second, existing methods employ a uniform, simplified approach for mating surfaces with different geometric features. The shape error of annular end faces exhibits a two-dimensional distribution within the circular region, its mathematical characteristic being an isotropic, diffuse deviation within the plane. The shape error of cylindrical mating surfaces shows a gradually varying trend along the axial direction and a periodic pattern along the circumference, its mathematical characteristic being a separable composite deviation in both the axial and circumferential directions. The fundamental difference in the spatial distribution patterns of these two types of surfaces necessitates different basis function systems for mathematical representation. Existing methods fail to differentiate modeling based on geometric feature types, resulting in low efficiency in deviation space description and potential omission of key morphological features.

[0008] Third, existing methods suffer from bottlenecks in computational efficiency regarding elastic contact deformation. Considering surface micro-morphology, the contact problem degenerates from ideal surface-to-surface contact into a collection of multiple discrete micro-protrusions. The contact area evolves asymptotically with the load, and the computational load increases quadratically with the number of mesh nodes. If traditional finite element methods are used to handle such problems, a single solution often takes tens of minutes to several hours, making it impossible to embed in Monte Carlo simulation loops requiring thousands of repetitive calculations. This forces existing methods to compromise between physical accuracy and computational feasibility.

[0009] Fourth, existing methods lack a quantitative decomposition of the contribution weights of each deviation source. Monte Carlo simulations ultimately provide the statistical distribution of assembly accuracy (mean and standard deviation, etc.), but do not further answer the questions truly needed for tolerance optimization decisions: which tolerance parameters contribute most to assembly accuracy fluctuations, which tolerances tightening yield the most significant improvement in accuracy, and which tolerances have reached saturation near their current values ​​and can be appropriately relaxed. This lack of quantitative sensitivity ranking information makes it difficult to effectively translate simulation results into process optimization decisions.

[0010] To address the aforementioned shortcomings, the purpose of this invention is to provide a Monte Carlo prediction method for assembly accuracy that considers microstructure and contact deformation. This method can establish parameterized models of non-ideal surface morphology that match the geometric characteristics of annular end faces and cylindrical mating surfaces; it can complete the calculation of elastic contact deformation considering surface microstructure in seconds during each Monte Carlo sampling, without needing to call the full finite element solution each time; it can quantitatively decompose the independent contribution ratio of each deviation source to the variation in assembly accuracy and output a sensitivity list sorted by contribution rate, providing a priority-based quantitative basis for tolerance optimization. Summary of the Invention

[0011] To address the shortcomings of existing technologies, the purpose of this invention is to provide a Monte Carlo prediction method for assembly accuracy that considers microstructure and contact deformation, enabling efficient statistical prediction and quantification of the contribution of deviation sources.

[0012] The technical solution adopted by this invention to solve its technical problem is:

[0013] A Monte Carlo method for predicting assembly accuracy that considers microstructure and contact deformation includes the following steps:

[0014] Step S1: Model simplification and finite element simulation. The assembly structure is geometrically cleaned and the model is simplified. Process features that do not affect the overall stiffness distribution are removed. A finite element model is constructed and boundary conditions and preloads are applied for solution. The contact pressure distribution of key mating surfaces under nominal working conditions is extracted.

[0015] Step S2: Model the non-ideal surface of the assembly interface, extract the geometric feature type and pressure distribution of the key mating surfaces based on the finite element model. The key mating surfaces include the annular end face and the cylindrical mating surface; establish the macroscopic pose deviation parameters and microscopic morphology deviation parameters respectively.

[0016] Step S3: Solve the contact deformation. For any combination of macroscopic pose deviation parameters and microscopic morphology deviation parameters, and in conjunction with the contact pressure distribution, solve the contact deformation of the key mating surfaces under the action of preload to obtain the contact state of the key mating surfaces.

[0017] Step S4: Axis pose fitting. Fit the spatial pose of the axis after assembly based on the contact state, and calculate the assembly accuracy index value from the axis spatial pose.

[0018] Step S5: Monte Carlo simulation. Repeat steps S2 to S4 for a total of N times to obtain N sets of assembly accuracy index values, forming a simulation sample dataset; where N is the preset total number of simulations.

[0019] Step S6: Result analysis. Perform statistical analysis on the simulation sample dataset, and output the statistical distribution characteristics of each assembly accuracy index and the sensitivity ranking of each macroscopic pose deviation parameter.

[0020] As a preferred embodiment, a further technical solution of the present invention is:

[0021] Preferably, the model simplification and finite element simulation in step S1 specifically include:

[0022] The assembly structure was geometrically cleaned and simplified, ignoring the thread structure of the bolts and simplifying them into smooth cylindrical structures. Process chamfers and fillets on the platform and inner sun ring that did not affect the overall stiffness distribution were removed. Mesh generation was performed using mesh mapping and projection methods, with variational meshes used appropriately to ensure the representation of circular hole features and the consistency of the overall mesh. The platform, inner ring, and intermediate shaft were all meshed using hexahedral structured meshes. Three contact pairs were defined: the intermediate shaft to the platform cylindrical surface, the intermediate shaft to the platform toroidal surface, and the intermediate shaft to the inner sun ring; all contact types were frictional. Material properties were set. Fixed constraints were applied to the bottom surfaces of the platform and inner sun ring, and a rated torque was applied to the circumferential threaded holes to simulate the bolt tightening process. The stress values ​​at each node on the cylindrical and toroidal surfaces were extracted as the contact pressure distribution.

[0023] Preferably, in step S2, the macroscopic pose deviation parameters are described using a small displacement screw parameter system. The small displacement screw parameters include three translation parameters and three rotation parameters. Each parameter is truncated within the tolerance zone and is randomly selected according to the probability distribution based on the statistical characteristics of the machining process. The parameters follow a normal distribution, and the standard deviation is calculated from the tolerance value marked on the drawing according to the 3σ principle.

[0024] Preferably, in step S2, the micro-morphological deviation parameters are modeled differently based on the geometric feature types of the key mating surfaces. For the annular end face, a Zernike polynomial is used to construct its shape error model. The Zernike polynomial is composed of the product of radial and angular functions. The shape deviation of any annular end face relative to the ideal plane is represented as a weighted superposition of several Zernike basis functions. The Zernike polynomial is truncated to the 16th order. For the cylindrical mating surface, a Legendre-Fourier polynomial is used to construct its shape error model. The Legendre-Fourier polynomial separates the cylindrical surface morphology into variables along the axial and circumferential directions. The axial component is described by the Legendre polynomial along the axial direction, and the circumferential component is described by the Fourier series along the circumferential direction. The actual radius deviation at any point on the cylindrical mating surface is represented as a weighted superposition of the combined product of the axial Legendre mode and the circumferential Fourier frequency. The Legendre polynomial is taken up to the 4th order, and the Fourier series is taken up to the 3rd order, which can characterize the axial bending, taper error, and radial elliptic and triangular errors of the cylindrical surface.

[0025] Preferably, the Zernike polynomial is defined in polar coordinates. On the unit disk below, the expression is:

[0026] ;

[0027] in, For normalized radial coordinates, , The actual distance from a point on the annular end face to the center of the circle. Where is the radius of the annular end face, 0 ≤ ≤1; Here, azimuth coordinates represent the angular position of a point on the annular end face in the circumferential direction, 0 ≤ 0. < ; The order of the radial polynomial controls the spatial frequency of the error. Let be the angular frequency parameter, satisfying | |≤ and -| | is an even number; It is a radial polynomial, and its expression is:

[0028] ;

[0029] in, The summation index takes values ​​of 0, 1, 2, ... -| | / 2

[0030] It is an angular function, and its expression is:

[0031] ;

[0032] The Legendre-Fourier polynomial in cylindrical coordinates is expressed as follows:

[0033] ;

[0034] in, Let be the angular coordinate, representing the angular position of a point on the mating surface of the cylinder in the circumferential direction, 0 ≤ 0. < ; The coordinates are axial, representing the position of a point on the mating surface of a cylinder along the axial direction. ∈[-1,1]; For non-ideal cylindrical surfaces in angular coordinates and axial coordinates The actual radius deviation at that location; For the first Legendre polynomials; and For the first Angular harmonic basis functions; and For joint expansion coefficients; The order index of the Fourier series. , The order index of Legendre polynomials , The maximum order of the preset Fourier series, The maximum order of the Legendre polynomial is preset;

[0035] The expression is:

[0036] ;

[0037] in, The summation index takes values ​​of 0, 1, 2, ..., [ ] indicates less than or equal to The largest integer.

[0038] Preferably, the contact deformation solution in step S3 is implemented using the conjugate gradient-fast Fourier transform method, specifically including: under the assumptions of elasticity, frictionlessness, and no adhesion, transforming the contact problem between two rough surfaces into the contact problem between a rigid smooth plane and an equivalent elastic rough surface, and the equivalent surface height function. From the formula Calculate, where, It is an equivalent rough surface height function. The surface morphology of the first contact body. The surface morphology of the second contact body is determined; the pressure distribution is iteratively searched using the conjugate gradient method. In each iteration, the spatial domain convolution is converted to the frequency domain dot product using the fast Fourier transform to accelerate the computation, and the contact pressure distribution obtained in step S1 is used as the initial value for the conjugate gradient method iteration; the equivalent elastic modulus is given by the formula... Calculate, where, For equivalent Young's modulus, , These are the Young's moduli of the two contacting bodies, respectively. , These are the Poisson's ratios of the two contacting bodies.

[0039] Preferably, the axis pose fitting in step S4 specifically involves: extracting the normal deformation of each discrete node on the key mating surface based on the contact state obtained in step S3; performing least-squares plane fitting on the node deformation in the contact area of ​​the annular end face to obtain the actual spatial pose of the annular end face; performing least-squares cylindrical surface fitting on the node deformation in the contact area of ​​the cylindrical mating surface to obtain the actual spatial axis of the cylindrical mating surface; using the actual spatial axis of the cylindrical mating surface as the spatial pose of the axis after assembly, calculating the positional deviation of the axis in the direction perpendicular to the axis, the tilt angle deviation in the mutually perpendicular projection planes, and the coaxiality deviation, and using the positional deviation, tilt angle deviation, and coaxiality deviation as assembly accuracy index values.

[0040] Preferably, the multi-process parallel computing architecture used in step S5 is as follows: Step S501: Pre-load geometric parameters, tolerance parameters, and the total number of simulations N. Based on the normal distribution characteristics of each macroscopic pose deviation parameter within the tolerance band, generate N sets of random sample data at once; Step S502: Distribute the N sets of random sample data to multiple independent computing processes. Each computing process executes the process from step S2 to step S4 once to obtain a set of assembly accuracy index values; Step S503: After all computing processes are completed, summarize the assembly accuracy index values ​​of each set to form a simulation sample dataset.

[0041] Preferably, the statistical distribution characteristics in step S6 include the mean. Standard deviation The KS test p-value is used to determine the normality of each assembly accuracy index. When the KS test p-value is greater than the preset significance level, the index is determined to follow a normal distribution. When the KS test p-value is less than or equal to the preset significance level, the index is determined to deviate from a normal distribution.

[0042] Preferably, in step S6, the sensitivity ranking is obtained by quantifying the independent contribution ratio of each macroscopic pose deviation parameter to the variation in assembly accuracy through variance contribution decomposition. Specifically, this includes: calculating the Pearson linear correlation coefficient between each macroscopic pose deviation parameter and the assembly accuracy index to reflect the degree of linear correlation between a single macroscopic pose deviation parameter and the assembly accuracy index; decomposing the total variance of the assembly accuracy index into the sum of the independent contribution term and the residual term of each macroscopic pose deviation parameter through variance decomposition, and calculating the independent contribution ratio of each macroscopic pose deviation parameter to the variation in the assembly accuracy index; generating a sensitivity ranking list in descending order of independent contribution ratio, the sensitivity ranking list including the independent contribution ratio value of each macroscopic pose deviation parameter.

[0043] The present invention, which adopts the above technical solution, has the following prominent features compared with the prior art:

[0044] (1) A parameterized model of micro-morphological deviation matching the geometric feature type was established. For the differences in geometric feature types of key mating surfaces, micro-morphological deviation parameters were constructed respectively. For the annular end face, a Zernike polynomial was used to construct a shape error model, utilizing its orthogonality and rotational invariance in the circular domain for morphological characterization. For the cylindrical mating surface, a Legendre-Fourier polynomial was used to construct a shape error model, separating the cylindrical surface morphology along the axial and circumferential directions. The axial component was described using a Legendre polynomial, and the circumferential component was described using a Fourier series, representing the periodic shape deviation along the circumference. This modeling method ensures that the mathematical structure of the micro-morphological deviation parameters matches the error intrinsics of the geometric features, facilitating the control of parameter dimensions while maintaining the accuracy of morphological description.

[0045] (2) A contact deformation solution mechanism based on contact pressure distribution and frequency domain transformation is formed. In the contact deformation solution, the contact pressure distribution obtained in step S1 is used as the initial value for iteration, and the conjugate gradient-fast Fourier transform method is used to accelerate the operation. By introducing the contact pressure distribution of the key mating surface extracted in step S1 under the nominal working condition into the iteration process of step S3, the convergence speed of the contact pressure distribution is effectively improved. At the same time, in each iteration, the fast Fourier transform is used to convert the spatial domain convolution into the frequency domain dot product, which reduces the computational complexity. This mechanism avoids running the complete finite element simulation repeatedly in each sampling of the Monte Carlo simulation, so that the single contact deformation solution can be completed in a short time, supporting the feasibility of repeating steps S2 to S4 in step S5.

[0046] (3) Sensitivity ranking output based on variance contribution decomposition is realized; In the result analysis, this application quantifies the independent contribution ratio of each macroscopic pose deviation parameter to the variation of assembly accuracy through variance contribution decomposition, and outputs the sensitivity ranking. Through the statistical analysis method in step S6, the total variance of the assembly accuracy index is decomposed into the sum of the independent contribution term and the residual term of each macroscopic pose deviation parameter, the independent contribution ratio of each macroscopic pose deviation parameter is calculated, and a sensitivity ranking list is generated according to the independent contribution ratio from large to small. This mechanism can identify macroscopic pose deviation parameters with high weights on the variation of assembly accuracy index, provide a decision basis for tolerance optimization, and help to prioritize the control of the tolerance range of key parameters in process design.

[0047] (4) A multi-process execution architecture for parallel Monte Carlo simulation is provided. In the parallel Monte Carlo simulation, this application adopts a multi-process parallel computing architecture to repeatedly execute steps S2 to S4. By pre-generating random sample data and distributing the samples to multiple independent computing processes, each computing process executes steps S2 to S4 once. Combined with the timeout monitoring mechanism within the computing process, the stability of the simulation process is ensured. This architecture effectively utilizes computing resources, enabling Monte Carlo simulations of a scale of thousands to be completed within an acceptable time. At the same time, the validity of the simulation sample dataset is ensured by skipping invalid samples. Attached Figure Description

[0048] Figure 1 This is a flowchart illustrating the overall process of the Monte Carlo prediction method for assembly accuracy in this embodiment of the invention.

[0049] Figure 2 This is a flowchart of the Monte Carlo simulation steps in an embodiment of the present invention;

[0050] Figure 3 This is a schematic diagram of the inertial platform assembly structure in an embodiment of the present invention;

[0051] Figure 4 This is a partially enlarged schematic diagram of a key assembly interface in an embodiment of the present invention;

[0052] Figure 5 This is a contact pressure distribution cloud map of a mating surface in an embodiment of the present invention;

[0053] Figure 6 This is a stress distribution diagram of the annular end face in an embodiment of the present invention;

[0054] Figure 7 This is a Zernike polynomial mode decomposition diagram of the microstructure of the annular end face in an embodiment of the present invention.

[0055] Figure 8 This is a schematic diagram of the Legendre-Fourier polynomial mode combination of the microstructure of the cylindrical mating surface in an embodiment of the present invention.

[0056] Figure 9 The axial position deviation in the embodiment of the present invention and Scatter plot;

[0057] Figure 10 This is a probability distribution histogram of assembly accuracy indicators in an embodiment of the present invention, where (a) and (b) represent the axis position deviations, respectively. and The normal distribution fitting curves, (c) and (d) are the axis tilt angle deviations, respectively. The distribution characteristic curve of the coaxiality deviation;

[0058] Figure 11 This is a heatmap of the Pearson linear correlation coefficients between various macroscopic pose deviation parameters and assembly accuracy indicators in this embodiment of the invention.

[0059] Figure 12 This is a bar chart showing the independent contribution ratio of each macroscopic pose deviation parameter to the coaxiality deviation variance in the embodiments of the present invention. Detailed Implementation

[0060] The present invention will be further illustrated below with reference to specific embodiments. The purpose of this illustration is solely to provide a better understanding of the invention. Therefore, the examples given do not limit the scope of protection of the present invention.

[0061] This invention provides a method for predicting the assembly accuracy of an inertial platform based on non-ideal surface morphology and Monte Carlo simulation. This method, by considering the microscopic geometric deviations of the contact surfaces and combining them with macroscopic assembly parameters, achieves high-precision prediction of accuracy indicators such as the axial position, tilt angle, and coaxiality after assembly. Figure 1 The diagram shown is the overall flowchart of the assembly accuracy Monte Carlo prediction method in this embodiment of the invention, illustrating the complete steps from model simplification, non-ideal surface modeling, contact deformation solving, axis pose fitting to Monte Carlo simulation and result analysis.

[0062] Furthermore, such as Figure 2 The flowchart shown is a Monte Carlo simulation process in an embodiment of the present invention, which specifically illustrates the logical relationship between random sample generation, sample allocation, single simulation process execution, and result summarization under a multi-process parallel computing architecture.

[0063] Step S1: Model simplification and finite element simulation. The assembly structure undergoes geometric cleanup and model simplification, removing process features that do not affect the overall stiffness distribution. A finite element model is constructed, and boundary conditions and preloads are applied for solution. The contact pressure distribution of key mating surfaces under nominal working conditions is extracted. For example... Figure 3 As shown, the inertial platform assembly structure mainly includes an intermediate shaft, a platform body, and an inner solar ring, such as... Figure 4 As shown, the geometric features of the two key contact areas—the intermediate axis and the cylindrical surface of the frustum, and the intermediate axis and the toroidal surface of the frustum—are further revealed through a magnified view.

[0064] Model simplification and finite element simulation specifically include:

[0065] The assembly structure underwent geometric cleanup and model simplification. To improve computational efficiency, the complex structure of the platform was simplified by ignoring the thread structure of the bolts and simplifying them into smooth cylindrical structures. Process chamfers and fillets on the platform and inner sun ring that do not affect the overall stiffness distribution were removed.

[0066] The mesh was generated using mesh mapping and projection methods, and variational meshing was appropriately applied to ensure the representation of the circular hole features and the consistency of the overall mesh. The platform, inner ring, and intermediate axis were all meshed using hexahedral structured meshes with Solid185 elements, totaling 153,904 meshes. Simultaneously, the mesh quality was checked to ensure that over 95% of the hexahedral elements had a Jacobian greater than 0.7, and the minimum Jacobian was no less than 0.47. In contact areas with large local stress gradients, the element aspect ratio was controlled within 5:1, and in transition regions, it did not exceed 8:1. The overall mesh quality met the solution requirements.

[0067] In the simulation parameter settings and model solution, the contact was first defined. The finite element model set three contact pairs: the contact between the intermediate shaft and the cylindrical surface of the platform, the contact between the intermediate shaft and the toroidal surface of the platform, and the contact between the intermediate shaft and the inner solar ring. All contact types were frictional contacts with a friction coefficient of 0.15. Each contact pair consisted of a contact surface (CONTA173) and a target surface (TARGE170).

[0068] The material properties are set as shown in Table 1.

[0069] Table 1 Material properties of the finite element model

[0070]

[0071] A fixed constraint is applied to the platform and the bottom surface of the inner solar ring. A torque of 2.2 N·m (the rated torque in this embodiment) is applied to the circumferential threaded hole to simulate the bolt tightening process. After the solution is completed, as follows... Figure 5 The color scale (right side) indicates the contact pressure magnitude (unit: MPa), visually presenting the stress concentration and distribution pattern on the mating surface; the stress values ​​of each node on the cylindrical mating surface and the annular end face are extracted as the contact pressure distribution, where the maximum stress on the cylindrical surface occurs at the location where the equivalent preload is applied, with a maximum value of 486.94 MPa; Figure 6 The stress distribution (unit: MPa) on the annular end face is displayed using a color scale (left side), with red circles marking key analysis areas. The maximum stress on the annular end face also corresponds to the location where the equivalent preload is applied, with a maximum value of 74.65 MPa. The stress values ​​at each node of the two key mating surfaces, the cylindrical surface and the annular end face, are used as load inputs for subsequent contact mechanics solutions.

[0072] Step S2: Model the non-ideal surface of the assembly interface. Extract the geometric feature type and pressure distribution of the key mating surfaces based on the finite element model. The key mating surfaces include the annular end face and the cylindrical mating surface. Establish macroscopic pose deviation parameters and microscopic morphology deviation parameters respectively.

[0073] (a) Macroscopic pose deviation parameters;

[0074] The macroscopic pose deviation parameters are described using the Small Displacement Torsor (SDT) parameter system. The Small Displacement Torsor parameters include three translation parameters and three rotation parameters, totaling six degrees of freedom deviation parameters, which are used to characterize the spatial pose deviation of each mating surface relative to its nominal design position.

[0075] Each parameter is truncated within the tolerance zone and randomly selected according to the probability distribution based on the statistical characteristics of the processing technology; each small displacement screw parameter follows a truncated normal distribution with a mean of zero, and a mean of 0 indicates no systematic deviation; the standard deviation of each parameter is determined by converting the tolerance value marked on the drawing according to the 3σ principle, that is, the relationship between the tolerance T and the standard deviation σ is T=6σ, σ=T / 6.

[0076] For annular end faces, when simultaneously subject to positioning tolerances With orientation tolerance When constrained, rotation parameters , The standard deviation is Translation parameters The standard deviation is ,in It is the outer diameter of the toroid.

[0077] For cylindrical mating surfaces, when simultaneously subject to positioning tolerances With orientation tolerance When constrained, rotation parameters , The standard deviation is Translation parameters , The standard deviation is ,in The height is the cylindrical surface.

[0078] The distribution settings of each parameter are shown in Table 2. All parameters are within the tolerance zone and follow a truncated normal distribution. The cutoff boundary is determined by the standard deviation and the tolerance range.

[0079] Table 2 SDT Parameter Probability Distribution Settings

[0080]

[0081] (ii) Microscopic morphology deviation parameters;

[0082] Microscopic morphology deviation parameters are modeled differently based on the geometric feature type of the key mating surfaces.

[0083] (1) Microscopic morphology modeling of the annular end face: The shape error model is constructed using Zernike polynomials, which are defined in polar coordinates. On the unit disk below, the expression is:

[0084] ;

[0085] in, For normalized radial coordinates, , The actual distance from a point on the annular end face to the center of the circle. Where is the radius of the annular end face, 0 ≤ ≤1; Here, azimuth coordinates represent the angular position of a point on the annular end face in the circumferential direction, 0 ≤ 0. < ; The order of the radial polynomial controls the spatial frequency of the error. Let be the angular frequency parameter, satisfying | |≤ and -| | is even to ensure rotational symmetry; It is a radial polynomial, and its expression is:

[0086] ;

[0087] in, The summation index takes values ​​of 0, 1, 2, ... -| | / 2

[0088] It is an angular function, and its expression is:

[0089] ;

[0090] Zernike polynomials possess orthogonality, completeness, and rotation invariance over the unit circle. Their orthogonality ensures that error modes of different orders are mathematically independent, which is crucial for error decomposition and coefficient solving; their rotation invariance allows for changes in the coordinate system when rotated by an angle. Afterwards, the function form remains unchanged.

[0091] Therefore, the Zernike polynomial can be expressed as:

[0092] ;

[0093] The above formula provides a schematic diagram of the first 16 terms of the Zernike polynomial basis functions, as shown below. Figure 7As shown, the shape deviation of an arbitrary annular end face relative to an ideal plane is represented as a weighted superposition of several Zernike basis functions, each corresponding to an independent error mode. The order of the basis functions determines the spatial resolution of the described morphology: lower-order terms correspond to global morphological deviations, while higher-order terms correspond to local texture details. In this embodiment, the Zernike polynomial is truncated to the 16th order, which can explain approximately 95% of the shape error variance of the annular end face, achieving a balance between description accuracy and computational efficiency.

[0094] (2) Microscopic morphology modeling of cylindrical mating surfaces: For cylindrical mating surfaces, the shape error model is constructed using Legendre-Fourier polynomials. Any cylindrical surface shape can be decomposed into a linear superposition of a set of Legendre-Fourier polynomials, where each polynomial corresponds to a specific error source in the machining process (such as clamping deformation, thermal drift, tool wear, etc.), which can systematically characterize the shape error of cylindrical surface features from the perspective of physical mechanism.

[0095] Legendre-Fourier polynomials separate the cylindrical surface morphology into variables along the axial and circumferential directions. The axial component is described by Legendre polynomials, representing the shape deviation along the axial direction, while the circumferential component is described by Fourier series, representing the periodic shape deviation along the circumference. The actual radius deviation at any point on the cylindrical mating surface is expressed as a weighted superposition of the combined product of the axial Legendre mode and the circumferential Fourier frequency.

[0096] Legendre polynomials are defined on intervals The q-th Legendre polynomial is denoted as... In the weight function When, the explicit expansion form is

[0097] ;

[0098] in, The summation index takes values ​​of 0, 1, 2, ..., [ ] indicates less than or equal to The largest integer.

[0099] Legendre polynomials possess important properties such as orthogonality and parity: different orders of Legendre polynomials satisfy orthogonality on the interval [-1,1], ensuring the independence between error components of each order; even-numbered terms are even functions and odd-numbered terms are odd functions, which can respectively characterize symmetrical errors (such as barrel shape and waist drum shape) and asymmetrical errors (such as cone shape and tilt) on the axial direction of the cylinder.

[0100] Legendre polynomials are suitable for expressing the axial errors of cylindrical features, while the radial (cross-sectional direction) errors of a cylinder generally exhibit periodicity. In engineering, the most ideal mathematical tool for describing such periodic distributions is the Fourier series. Therefore, this application uses a combination of Legendre polynomials and Fourier series to model cylindrical errors, with the basic form as follows:

[0101] ;

[0102] In the formula, A periodic function representing the circumferential direction (such as the angular distribution of roundness / ripple error). Angular coordinates, This is a constant term (average component). , Let N be the Fourier series of the p-th cosine and sine terms, where p is the harmonic order and N is the order of the Fourier series.

[0103] In cylindrical error modeling, Fourier series is used to describe radial ( Periodic errors in direction, such as roundness error, waviness, and polygonal distortion. By combining Legendre polynomials with Fourier series, a Legendre-Fourier polynomial (e.g., Figure 8 As shown in the figure, it is used to simultaneously describe the combined axial and radial errors of a cylindrical surface feature. Its expression in the cylindrical coordinate system is:

[0104] ;

[0105] in, Let be the angular coordinate, representing the angular position of a point on the mating surface of the cylinder in the circumferential direction, 0 ≤ 0. < ; The coordinates are axial, representing the position of a point on the mating surface of a cylinder along the axial direction. ∈[-1,1]; For non-ideal cylindrical surfaces in angular coordinates and axial coordinates The actual radius deviation at that location; For the first Legendre polynomials; and For the first Angular harmonic basis functions; and For joint expansion coefficients; The order index of the Fourier series. , The order index of Legendre polynomials , The maximum order of the preset Fourier series, This is the maximum order of the Legendre polynomial as preset.

[0106] By appropriately selecting the truncation order and solving the coefficient matrix, a precise mathematical description of the characteristic shape error of the cylinder can be obtained. In a specific embodiment, the Legendre polynomial is taken up to the 4th order and the Fourier series is taken up to the 3rd order, which can effectively characterize the axial bending and taper errors of the cylinder and the radial elliptic and triangular errors.

[0107] Step S3: Solve for contact deformation. For any combination of macroscopic pose deviation parameters and microscopic morphology deviation parameters, and in conjunction with the contact pressure distribution, solve for the contact deformation of the key mating surfaces under the action of preload, and obtain the contact state of the key mating surfaces.

[0108] The contact deformation is solved using the conjugate gradient-fast Fourier transform (CG-FFT) method. Specifically, under the assumptions of elasticity, frictionlessness, and no adhesion, the contact problem between two rough surfaces is transformed into a contact problem between a rigid smooth plane and an equivalent elastic rough surface. The equivalent surface height function... From the formula Calculate, where, It is an equivalent rough surface height function. The surface morphology of the first contact body. This represents the surface morphology of the second contact body.

[0109] When the original surface is expanded using a unified basis function, the equivalent surface can be represented as:

[0110] ;

[0111] in, and The average height of the two surfaces, , For expansion coefficients, The expression uses the basis functions. This expression transforms the complex three-dimensional contact problem of two rough surfaces into a single-surface contact model between an equivalent elastic body and a rigid plane, laying the theoretical foundation for subsequent numerical calculations. Under external loads, the normal displacements of each discrete node on the equivalent elastic surface are given by the following equation through convolution;

[0112] ;

[0113] in, For nodal displacement, and Discrete meshes in and Number of nodes in the direction, For nodes Contact pressure at the point, The influence coefficient function characterizes the elastic displacement response caused by a unit force, and its specific expression is given by the following equation:

[0114] ;

[0115] ;

[0116] and These represent the discrete mesh in... and Discrete spacing in direction, and For relative node indexes, Equivalent Young's modulus;

[0117] Equivalent Young's modulus The elastic modulus and Poisson's ratio of the two contacting materials are synthesized according to the following formula:

[0118] ;

[0119] in, For equivalent Young's modulus, , These are the Young's moduli of the two contacting bodies, respectively. , These are the Poisson's ratios of the two contacting bodies.

[0120] Under non-adhesive contact conditions, the system must satisfy the following contact constraints:

[0121] Non-penetrating constraint: The node contact gap is greater than or equal to zero, meaning that the surfaces cannot penetrate each other after contact;

[0122] Non-negative pressure constraint: The contact pressure is always compressive stress, meaning the node cannot withstand tensile stress;

[0123] Complementary conditions: The nodes are either in contact (zero gap, positive pressure) or in non-contact (positive gap, zero pressure).

[0124] Based on the principle of minimum potential energy, this contact problem can be transformed into the following constrained optimization problem:

[0125] ;

[0126] in, For contact pressure vector, M is the total number of discrete nodes. Let the initial gap vector be... , for Influence coefficient matrix.

[0127] To efficiently solve this optimization problem, the conjugate gradient method is used for iterative solution, with the contact pressure distribution obtained in step S1 serving as the initial value for the conjugate gradient method iteration. Since this initial value is already within a physically reasonable neighborhood, the number of iterations required for convergence is significantly reduced.

[0128] In each iteration, the pressure distribution estimate is updated along the negative gradient direction, and the step size that minimizes the total residual energy is searched along this direction. Then, the pressure distribution is updated, and the convergence condition is checked. The convergence criterion is that the global maximum value of the relative change in pressure at each node between two consecutive iterations is less than a preset threshold. In this embodiment, the preset threshold is set to 10. -5 In practical applications, the threshold can be 10. -4 Up to 10 -6 The threshold is selected within a certain range. The smaller the threshold, the higher the calculation accuracy but the more iterations are required. The larger the threshold, the faster the calculation speed but the lower the accuracy. Those skilled in the art can determine the threshold within this range according to specific accuracy requirements.

[0129] In each iteration, the elastic displacement response is calculated using the convolution of the influence coefficient matrix and the pressure distribution. This convolution operation is accelerated in the frequency domain using a Fast Fourier Transform (FFT): first, the pressure distribution and the influence coefficient matrix are each subjected to a two-dimensional Fourier transform; then, the corresponding frequency components are multiplied point-by-point in the frequency domain; finally, the product is subjected to an inverse Fourier transform to obtain the displacement distribution. Since frequency domain multiplication is equivalent to spatial domain convolution, and the arithmetic complexity of frequency domain operations is the number of nodes multiplied by the logarithm of the number of nodes, the computation time is reduced from hours to seconds when the number of grid nodes reaches hundreds of thousands.

[0130] By solving the above CG-FFT iteratively, the contact pressure distribution and contact gap of each node of the key mating surface are obtained, that is, the contact state of the key mating surface.

[0131] Step S4: Axis pose fitting. Fit the spatial pose of the assembled axis according to the contact state, and calculate the assembly accuracy index value from the axis spatial pose.

[0132] In step S4, the axis pose fitting is specifically as follows: Based on the contact state obtained in step S3, the normal deformation of each discrete node on the key mating surface is extracted; the node deformation in the contact area of ​​the annular end face is fitted with least squares plane to obtain the actual spatial pose of the annular end face; the node deformation in the contact area of ​​the cylindrical mating surface is fitted with least squares cylindrical surface to obtain the actual spatial axis of the cylindrical mating surface; using the actual spatial axis of the cylindrical mating surface as the spatial pose of the axis after assembly, the positional deviation of the axis in the direction perpendicular to the axis, the tilt angle deviation in the mutually perpendicular projection plane, and the coaxiality deviation are calculated, and the positional deviation, tilt angle deviation, and coaxiality deviation are used as assembly accuracy index values.

[0133] Step S5: Monte Carlo simulation. Repeat steps S2 to S4 for a total of N times to obtain N sets of assembly accuracy index values, forming a simulation sample dataset; where N is the preset total number of simulations.

[0134] The statistical accuracy of Monte Carlo simulations is closely related to the sample size. According to the central limit theorem, the standard error of the mean estimate is... ,in The standard deviation of the estimated assembly accuracy index. The sample size is given by the formula [formula missing]. To achieve the preset accuracy, the required sample size is determined by [formula missing]. Confirmed, among which The quantiles of the standard normal distribution To account for the allowable estimation error, in this embodiment, considering both computational efficiency and accuracy, the total number of simulations N is set to 1000.

[0135] Each virtual assembly test is independent of the others and is executed using a multi-process parallel computing architecture, specifically including the following sub-steps:

[0136] Step S501: Pre-load geometric parameters, tolerance parameters, and the total number of simulations N. Based on the truncated normal distribution characteristics of each macroscopic pose deviation parameter within the tolerance band, generate N sets of random sample data at once. Each set of random sample data contains all the value combinations of macroscopic pose deviation parameters from step S2. Each parameter is independently extracted within its tolerance band according to a truncated normal distribution to avoid introducing bias due to repeated random number seeds or sequence correlation in a parallel environment. Among them, the macroscopic pose deviation parameters are independently extracted within their tolerance band according to a truncated normal distribution; the microscopic morphology deviation parameters include the expansion coefficients of Zernike polynomials and the joint expansion coefficients of Legendre-Fourier polynomials. Each expansion coefficient is randomly selected according to a normal or uniform distribution based on the statistical characteristics of the processed surface morphology.

[0137] Step S502: Distribute N sets of random sample data into multiple independent computation processes, each processing one or more samples independently. The computation process for each sample includes: extracting the macroscopic pose deviation parameters and microscopic shape deviation parameters corresponding to the current sample from the pre-generated random sample matrix; constructing the 3D shape of the mating surface according to the non-ideal surface modeling method in step S2; inputting the constructed non-ideal surface shape into the CG-FFT contact solver in step S3 for contact deformation calculation; and obtaining the assembly accuracy index value through axis pose fitting in step S4. Each computation process executes a complete process from step S2 to step S4 once, obtaining a set of assembly accuracy index values.

[0138] Step S503: To ensure computational robustness and efficiency, a timeout monitoring and retry mechanism is implemented within each computation process. If a single CG-FFT calculation fails to converge within the set time, the calculation for that sample is automatically restarted until convergence is achieved or the maximum number of retries is reached, thus preventing individual abnormal samples from blocking the overall parallel process. The computation results of each process are sent back to the main process for unified recording. After all computation processes are completed, the effective assembly accuracy index values ​​of each group are summarized, and invalid sample data from failed timeout retries are removed to form a simulation sample dataset.

[0139] Using the parallel Monte Carlo simulation framework described above, simulations of a thousand scale can be completed within an acceptable time, providing sufficient sample data for the statistical analysis in step S6.

[0140] Step S6: Result analysis. Perform statistical analysis on the simulation sample dataset, and output the statistical distribution characteristics of each assembly accuracy index and the sensitivity ranking of each macroscopic pose deviation parameter.

[0141] (I) Statistical Distribution Characteristics Analysis

[0142] Statistical analysis is performed on the simulation sample dataset obtained in step S5 to extract the probability distribution characteristics of each assembly accuracy index. The assembly accuracy index includes the axial position deviation in two orthogonal directions. and axial tilt deviation And coaxiality deviation. For example... Figure 9 As shown, it gives and The scatter distribution shows that most sample points are concentrated in the range of [-20μm, 20μm], and the scatter points are symmetrically distributed around the origin, without showing obvious directional systematic deviation. This indicates that the assembly error amplitudes in the x and y directions are comparable and have independent random fluctuation characteristics.

[0143] Probability and statistical distribution analysis was performed on the four deviation parameters mentioned above, and the results are as follows: Figure 10As shown, regarding the deviation of the axis position and Perform normal distribution fitting and plot the probability distribution histogram. Figure 10 In (a) and (b), the histograms and corresponding fitted curves match well, both exhibiting a typical bell-shaped symmetrical distribution, consistent with the theoretical expectation that random errors follow a normal distribution; simultaneously, the axial tilt deviation is analyzed. Distribution characteristics of coaxiality deviation ( Figure 10 (c) and (d) show a certain degree of skewness, indicating that they may not follow a simple normal distribution and require further testing.

[0144] Calculate the statistical characteristic values ​​of each assembly accuracy index, including the mean. Standard deviation And the KS test, which can describe the distribution pattern.

[0145] The normality of the four parameters was determined by the KS test, and the p-values ​​of the KS test are shown in Table 3.

[0146] Table 3 Statistical results of axis pose deviation

[0147]

[0148] Axis position deviation and The p-values ​​of all values ​​were greater than the significance level, so the null hypothesis that the distribution follows a normal distribution is accepted; while The p-value of the coaxiality deviation is lower than the significance threshold, indicating that its distribution deviates significantly from normality and is relatively complex.

[0149] Table 3 summarizes the mean, standard deviation, maximum and minimum values ​​of each assembly accuracy index. The mean value is -0.040 μm. The mean value is 0.000 μm, and the absolute values ​​are all much smaller than the corresponding standard deviations (8.400 μm and 8.440 μm), indicating that the assembly deviations are generally concentrated near the zero position and there is no significant systematic shift. The statistical values ​​for coaxiality deviation reflect a wider range of dispersion and a possible asymmetric distribution, which corresponds to the test results showing that it does not conform to a normal distribution; axis angle deviation The coaxiality deviation exhibits a significantly skewed distribution (non-normal), with KS test p-values ​​all less than 0.01. Preliminary analysis suggests the cause is that coaxiality deviation is defined as a non-linear combination of positional deviations in the x and y directions. Even... and If a particle follows a normal distribution with a mean of zero, the square root of its sum of squares will naturally follow a Rayleigh distribution (coaxiality = ...). (), exhibiting a skewed characteristic starting from 0 and having a long right tail. Similarly, The positive skewness may originate from the vector synthesis effect of the angular deviation, as well as the non-negative constraint of the cylindrical orientation tolerance zone on the rotation angle.

[0150] (II) Parameter Correlation Analysis

[0151] Calculate the Pearson linear correlation coefficient between each macroscopic pose deviation parameter and each assembly accuracy index. The closer the absolute value of the correlation coefficient is to 1, the stronger the linear correlation between the two variables; the closer it is to 0, the weaker the linear correlation. A significant correlation is determined by a correlation coefficient absolute value greater than 0.5. When |r|>0.5, the macroscopic pose deviation parameter and the assembly accuracy index are considered to have a significant correlation; when |r|≤0.5, the linear correlation is considered not significant. Through this screening, the set of parameters that play a dominant role in assembly accuracy can be quickly identified, and key parameters with significant correlation to assembly accuracy can be selected.

[0152] like Figure 11 As shown, the thermodynamic relationship between the small displacement screw parameter and the axial position deviation, tilt angle deviation, and coaxiality deviation is illustrated, reflecting the linear correlation between each tolerance element and assembly deviation. It can be seen that the axial position deviation... The translation parameter u2 of the cylinder along the x-direction is positively correlated with the height. The translation parameter v2 of the cylinder along the y-direction is also positively correlated with the height. Statistical results show that the translations u2 and v2 of the cylinder are related to the deviation of the axis position ( , The correlation with the axis inclination angle is the strongest, indicating that its variation has a direct impact on the alignment accuracy of the axis. All rotational parameters (α1, β1, α2, β2) showed a certain degree of correlation, indicating that small rotational changes of the cylindrical surface around the x and y axes would collectively affect the axial direction after assembly, and deviation of any single rotational element could be reflected in the angular offset of the axis. Furthermore, the axial translation parameter w1 of the toroidal surface had a significant impact on the axial position of the assembly, but its correlation coefficient with coaxiality deviation was low, indicating that simple axial movement had a relatively small impact on the coaxiality deviation, which was still mainly determined by lateral translation and angular deviation.

[0153] (III) Sensitivity Ranking

[0154] The contribution of the small displacement screw (SDT) parameter to the coaxiality variance is as follows: Figure 12The influence of different tolerance elements on assembly output variation was quantified by variance decomposition. It can be seen that the translation parameter v2 along the y-direction of the cylinder has the highest contribution rate, approximately 66%, indicating that under the current tolerance setting, this parameter is the most significant source of error causing coaxiality fluctuations. This is followed by the translation parameter u2 along the x-direction (approximately 10%), the rotation parameter α2 about the x-axis (approximately 10%), and the rotation parameter β2 about the y-axis (approximately 8%), which together contribute approximately 28%, constituting secondary influencing factors. In contrast, the total contribution rate of the toroidal parameters α1, β1, and w1 is only about 5%, indicating that the geometrical variation of the toroid has a very weak impact on coaxiality. However, w1 mainly affects the axial position of the assembly and has a relatively small impact on coaxiality. Therefore, the positioning accuracy of the cylindrical surface is a key factor in determining the coaxiality of the assembly. The translation in the v2 direction is particularly important, as even a small fluctuation can cause a significant change in coaxiality. In tolerance design, the positioning tolerance of the cylindrical surface should be optimized first, and the variation range of u2 and v2 should be strictly constrained to achieve precise control of the assembly coaxiality.

[0155] In summary, the parameterized model of non-ideal mating surfaces established in this invention, based on small displacement screw parameters and orthogonal polynomials, describes macroscopic pose deviations through small displacement screw parameters and constructs microscopic morphology deviation parameters for the annular end face and cylindrical mating surfaces using Zernike polynomials and Legendre-Fourier polynomials, respectively, effectively balancing the systematic deviations of assembly poses and the random fluctuations of surface morphology. The constructed parallel solution framework for assembly deviations with toroidal-cylindrical dual-interface coupling introduces contact pressure distribution as the initial value for iteration and uses the conjugate gradient-fast Fourier transform method to solve for contact deformation, significantly reducing the time consumption of a single simulation and supporting the feasibility of Monte Carlo simulations on a scale of thousands. The sensitivity ranking results based on variance contribution decomposition output show that the translation parameter along the y-axis of the cylindrical mating surface has the highest independent contribution ratio to the coaxiality deviation variance, providing a quantitative basis for prioritizing the control of key parameters and rationally allocating manufacturing resources in tolerance optimization.

[0156] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of the invention. All equivalent changes made based on the description and drawings of the present invention are included within the scope of the present invention.

Claims

1. A Monte Carlo method for predicting assembly accuracy considering microstructure and contact deformation, characterized in that, Includes the following steps: Step S1: Model simplification and finite element simulation. The assembly structure is geometrically cleaned and the model is simplified. Process features that do not affect the overall stiffness distribution are removed. A finite element model is constructed and boundary conditions and preloads are applied for solution. The contact pressure distribution of key mating surfaces under nominal working conditions is extracted. Step S2: Model the non-ideal surface of the assembly interface, extract the geometric feature type and pressure distribution of the key mating surfaces based on the finite element model. The key mating surfaces include the annular end face and the cylindrical mating surface; establish the macroscopic pose deviation parameters and microscopic morphology deviation parameters respectively. Step S3: Solve the contact deformation. For any combination of macroscopic pose deviation parameters and microscopic morphology deviation parameters, and in conjunction with the contact pressure distribution, solve the contact deformation of the key mating surfaces under the action of preload to obtain the contact state of the key mating surfaces. Step S4: Axis pose fitting. Fit the spatial pose of the axis after assembly based on the contact state, and calculate the assembly accuracy index value from the axis spatial pose. Step S5: Monte Carlo simulation. Repeat steps S2 to S4 for a total of N times to obtain N sets of assembly accuracy index values, forming a simulation sample dataset; where N is the preset total number of simulations. Step S6: Result analysis. Perform statistical analysis on the simulation sample dataset, and output the statistical distribution characteristics of each assembly accuracy index and the sensitivity ranking of each macroscopic pose deviation parameter.

2. The Monte Carlo method for predicting assembly accuracy considering microstructure and contact deformation according to claim 1, characterized in that: Step S1, specifically the model simplification and finite element simulation, includes: The assembly structure is geometrically cleaned and the model is simplified. The thread structure of the bolts is ignored, and the bolts are simplified into smooth cylindrical structures. The process chamfers and fillets on the platform and inner sun ring that do not affect the overall stiffness distribution are removed. The mesh was generated using mesh mapping and projection methods, and variational meshes were used appropriately to ensure the representation of the circular hole features and the coherence of the overall mesh. The platform, inner ring, and intermediate shaft were all meshed using hexahedral structured meshes. Three contact pairs were defined: the intermediate shaft and the cylindrical surface of the platform, the intermediate shaft and the toroidal surface of the platform, and the intermediate shaft and the inner sun ring. All contact types were frictional contacts. Material properties were set. Fixed constraints were applied to the bottom surfaces of the platform and the inner sun ring. A rated torque was applied to the circumferential threaded hole to simulate the bolt tightening process. The stress values ​​of each node on the cylindrical and toroidal surfaces were extracted as the contact pressure distribution.

3. The Monte Carlo prediction method for assembly accuracy considering microstructure and contact deformation according to claim 1, characterized in that: In step S2, the macroscopic pose deviation parameters are described using a small displacement screw parameter system. The small displacement screw parameters include three translation parameters and three rotation parameters. Each parameter is truncated within the tolerance zone and is randomly selected according to the probability distribution based on the statistical characteristics of the machining process. The parameters follow a normal distribution, and the standard deviation is calculated from the tolerance value marked on the drawing according to the 3σ principle.

4. The Monte Carlo prediction method for assembly accuracy considering microstructure and contact deformation according to claim 1, characterized in that: In step S2, the micro-morphological deviation parameters are modeled differently based on the geometric characteristics of the key mating surfaces. For the annular end face, a Zernike polynomial is used to construct its shape error model. The Zernike polynomial is composed of the product of radial and angular functions. The shape deviation of any annular end face relative to the ideal plane is represented as a weighted superposition of several Zernike basis functions. The Zernike polynomial is truncated to the 16th order. For the cylindrical mating surface, a Legendre-Fourier polynomial is used to construct its shape error model. The Legendre-Fourier polynomial separates the cylindrical surface morphology into variables along the axial and circumferential directions. The axial component is described by the Legendre polynomial along the axial direction, and the circumferential component is described by the Fourier series along the circumferential direction. The actual radius deviation at any point on the cylindrical mating surface is represented as a weighted superposition of the combined product of the axial Legendre mode and the circumferential Fourier frequency. The Legendre polynomial is taken up to the 4th order, and the Fourier series is taken up to the 3rd order, which can characterize the axial bending, taper error, and radial elliptic and triangular errors of the cylindrical surface.

5. The Monte Carlo method for predicting assembly accuracy considering microstructure and contact deformation according to claim 4, characterized in that: Zernike polynomials are defined in polar coordinates On the unit disk below, the expression is: ； in, For normalized radial coordinates, , The actual distance from a point on the annular end face to the center of the circle. Where is the radius of the annular end face, 0 ≤ ≤1; Here, azimuth coordinates represent the angular position of a point on the annular end face in the circumferential direction, 0 ≤ 0. < ; The order of the radial polynomial controls the spatial frequency of the error. Let be the angular frequency parameter, satisfying | |≤ and -| | is an even number; It is a radial polynomial, and its expression is: ； in, The summation index takes values ​​of 0, 1, 2, ... -| | / 2 It is an angular function, and its expression is: ； The Legendre-Fourier polynomial in cylindrical coordinates is expressed as follows: ； in, Let be the angular coordinate, representing the angular position of a point on the mating surface of the cylinder in the circumferential direction, 0 ≤ 0. < ; The coordinates are axial, representing the position of a point on the mating surface of a cylinder along the axial direction. ∈[-1,1]; For non-ideal cylindrical surfaces in angular coordinates and axial coordinates The actual radius deviation at that location; For the first Legendre polynomials; and For the first Angular harmonic basis functions; and For joint expansion coefficients; The order index of the Fourier series. , The order index of Legendre polynomials , The maximum order of the preset Fourier series, The maximum order of the Legendre polynomial is preset; The expression is: ； in, The summation index takes values ​​of 0, 1, 2, ..., [ ] indicates less than or equal to The largest integer.

6. The Monte Carlo method for predicting assembly accuracy considering microstructure and contact deformation according to claim 1, characterized in that, In step S3, the contact deformation is solved using the conjugate gradient-fast Fourier transform method, which specifically includes: Under the assumptions of elasticity, frictionlessness, and no adhesion, the contact problem between two rough surfaces is transformed into the contact problem between a rigid smooth plane and an equivalent elastic rough surface, where the equivalent surface height function is... From the formula Calculate, where, It is an equivalent rough surface height function. The surface morphology of the first contact body. The surface morphology of the second contact body; The pressure distribution is iteratively searched using the conjugate gradient method. In each iteration, the spatial domain convolution is converted into the frequency domain dot product using the fast Fourier transform to accelerate the operation. The contact pressure distribution obtained in step S1 is used as the initial value for the conjugate gradient method. The equivalent elastic modulus is given by the formula, Calculate, where, For equivalent Young's modulus, , These are the Young's moduli of the two contacting bodies, respectively. , These are the Poisson's ratios of the two contacting bodies.

7. The Monte Carlo method for predicting assembly accuracy considering microstructure and contact deformation according to claim 1, characterized in that, In step S4, the axis pose fitting is specifically as follows: based on the contact state obtained in step S3, the normal deformation of each discrete node on the key mating surface is extracted; the node deformation in the contact area of ​​the annular end face is fitted with least squares plane to obtain the actual spatial pose of the annular end face; the node deformation in the contact area of ​​the cylindrical mating surface is fitted with least squares cylindrical surface to obtain the actual spatial axis of the cylindrical mating surface. Using the actual spatial axis of the cylindrical mating surface as the spatial pose of the axis after assembly, the positional deviation of the axis in the direction perpendicular to the axis, the tilt angle deviation in the mutually perpendicular projection plane, and the coaxiality deviation are calculated. The positional deviation, tilt angle deviation, and coaxiality deviation are used as assembly accuracy index values.

8. The Monte Carlo prediction method for assembly accuracy considering microstructure and contact deformation according to claim 1, characterized in that, The multi-process parallel computing architecture used in step S5 is specifically as follows: Step S501: Pre-load geometric parameters, tolerance parameters, and total number of simulations N. Based on the normal distribution characteristics of each macroscopic pose deviation parameter within the tolerance band, generate N sets of random sample data at once. Step S502: Distribute N sets of random sample data into multiple independent calculation processes. Each calculation process executes the process from step S2 to step S4 once to obtain a set of assembly accuracy index values. Step S503: After all calculation processes are completed, summarize the assembly accuracy index values ​​of each group to form a simulation sample dataset.

9. The Monte Carlo method for predicting assembly accuracy considering microstructure and contact deformation according to claim 1, characterized in that, The statistical distribution characteristics in step S6 include the mean. Standard deviation The KS test p-value is used to determine the normality of each assembly accuracy index. When the KS test p-value is greater than the preset significance level, the index is determined to follow a normal distribution. When the KS test p-value is less than or equal to the preset significance level, the index is determined to deviate from a normal distribution.

10. The Monte Carlo method for predicting assembly accuracy considering microstructure and contact deformation according to claim 9, characterized in that, In step S6, the sensitivity ranking is obtained by quantifying the independent contribution ratio of each macroscopic pose deviation parameter to the variation in assembly accuracy through variance contribution decomposition, specifically including: The Pearson linear correlation coefficient between each macroscopic pose deviation parameter and the assembly accuracy index is calculated to reflect the degree of linear correlation between a single macroscopic pose deviation parameter and the assembly accuracy index. The total variance of the assembly accuracy index is decomposed into the sum of the independent contribution terms and residual terms of each macroscopic pose deviation parameter by variance decomposition, and the independent contribution ratio of each macroscopic pose deviation parameter to the variation of the assembly accuracy index is calculated. A sensitivity ranking list is generated in descending order of independent contribution ratio. The sensitivity ranking list includes the independent contribution ratio values ​​of each macroscopic pose deviation parameter.

Images