Fatigue life prediction method and system for high-strength bolt connection

Abstract

The invention relates to the technical field of bolt connection fatigue life prediction, in particular to a high-strength bolt connection fatigue life prediction method and system, and the method comprises the steps: obtaining a bolt entity morphology discrete point cloud, and generating a continuous geometric topology function through non-uniform rational B-spline interpolation reconstruction; discretizing the standard thread parameters to construct a standard finite element model, calculating a normal spatial position difference value based on a continuous geometric topology function, and packaging to generate a geometric defect mapping tensor; according to the converted node coordinates, constructing a real morphology finite element model; solving the nonlinear equilibrium equation by utilizing a penalty function method, and generating a contact stiffness characteristic matrix; calculating an axial load distribution sequence and executing gradient convolution operation, and extracting a local field variable peak value; and calculating an accumulated damage index in combination with a damage evolution equation and outputting a predicted life. By mapping the entity discrete point cloud into a finite element geometric defect tensor and reconstructing a real morphology model, the accurate prediction of the fatigue life of the bolt considering the micromachining deviation influence is realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of bolt connection fatigue life prediction, in particular to a high-strength bolt connection fatigue life prediction method and system. BACKGROUND

[0002] In the numerical simulation analysis of large bridge cable end anchoring systems, the existing calculation method usually uses a standard helical surface generated based on an analytical geometry equation for meshing, assuming that the node coordinates on the thread surface strictly follow the ideal mathematical trajectory. However, in the actual cable end anchoring scene, a large amount of discrete point cloud data on the bolt surface can be obtained through high-precision optical screening equipment.

[0003] However, the existing simulation algorithm cannot directly utilize these discrete point cloud data. The current mesh generation algorithm usually smooths or regularizes the input discrete data to adapt to the finite element topology structure, and this process loses the high-frequency local geometric deviation information (such as the small cumulative step of pitch or the random disturbance of the flank angle). Due to the lack of a calculation model that can map the discrete detection data to the finite element boundary conditions without loss, the solver cannot calculate the local singularity of the contact stiffness matrix caused by the microscopic geometric topology difference, thereby limiting the prediction accuracy of the fatigue life to the idealized assumption.

[0004] The existing technology cannot establish a direct mapping algorithm between the discrete appearance detection data and the continuous finite element mesh topology, resulting in the inability to preserve the microscopic geometric defect characteristics in numerical calculation, and thus it is difficult to accurately solve the problem of local non-uniform distribution of contact stiffness caused by processing deviation.

[0005] Therefore, the present application provides a high-strength bolt connection fatigue life prediction method and system. SUMMARY

[0006] The purpose of the present application is to provide a high-strength bolt connection fatigue life prediction method and system, which maps the entity discrete point cloud to the finite element geometric defect tensor and reconstructs the real topography model, realizing the accurate prediction of bolt fatigue life considering the influence of microscopic processing deviation.

[0007] To achieve the above purpose, the present application provides the following technical solutions: A high-strength bolt connection fatigue life prediction method, comprising: obtaining a discrete point cloud coordinate sequence representing the entity topography of the target batch of high-strength bolts; reconstructing the discrete point cloud coordinate sequence by a non-uniform rational B-spline interpolation method to generate a continuous geometric topology function; The preset standard thread geometry design parameters are discretized by a meshing method to construct a standard finite element mesh model; based on the continuous geometry topology function and the standard finite element mesh model, the normal space position difference is calculated and matrixed to generate a geometry defect mapping tensor; Based on the geometry defect mapping tensor, the node coordinate transformation is performed on the standard finite element mesh model to construct a real topography finite element model; based on the real topography finite element model, the penalty function method is used to construct a nonlinear numerical balance equation set, and the normal interaction value between each discrete contact node pair is iteratively solved to generate a contact stiffness feature matrix; Based on the contact stiffness feature matrix, the load distribution feature sequence distributed along the axial direction is calculated, and the gradient convolution operation is performed on the load distribution feature sequence to calculate the local field variable peak value; Based on the local field variable peak value, the cumulative damage index is calculated by combining the damage evolution equation, and the predicted life of the target batch of high-strength bolts is output.

[0008] Preferably, the process of obtaining the discrete point cloud coordinate sequence representing the real topography of the target batch of high-strength bolts comprises: The original three-dimensional scattered point cloud data output by the optical screening device is read, and a three-dimensional KD tree spatial index structure is constructed; each data point in the original three-dimensional scattered point cloud data is traversed, the K neighborhood points of the data point are searched using the three-dimensional KD tree spatial index structure, and the average Euclidean distance of the data point to all K neighborhood points is calculated; the mean and standard deviation of the average Euclidean distances of all data points are calculated, and the outlier noise points with an average Euclidean distance greater than the mean plus twice the standard deviation are removed; the remaining data points after removing the noise points are mapped from the Cartesian coordinate system to the cylindrical coordinate system, and the data points are rearranged in order according to the increasing phase angle value of the spiral line to generate an ordered discrete point cloud coordinate sequence.

[0009] Preferably, the generation process of the continuous geometry topology function comprises: The parameter node vector of each data point in the discrete point cloud coordinate sequence is calculated using the cumulative chord length parameterization method; a linear approximation equation set based on the least squares method is constructed, the discrete point cloud coordinate sequence is substituted into the linear approximation equation set as the target observation value, and the control vertex coordinates that minimize the equation set residual sum of squares are solved; the control vertex coordinates, the parameter node vector and the preset basis function order are used to construct a parameterized continuous geometry topology function, and the continuous geometry topology function is an analytical expression describing the geometry trajectory of the thread surface.

[0010] Preferably, the generation process of the geometry defect mapping tensor comprises: The standard thread geometry design parameter is read to generate an ideal spiral surface, and a hexahedral element is used to discretize the ideal spiral surface in space to construct a standard finite element grid model containing a node coordinate matrix and an element connection relationship; the discrete nodes on the contact surface in the standard finite element grid model are traversed, and a ray equation is constructed along the unit normal vector of each discrete node; the ray equation and the continuous geometric topology function are solved by a numerical root algorithm to obtain the spatial intersection coordinates of the analytical surface defined by the ray equation and the continuous geometric topology function; the Euclidean distance vector between the spatial intersection coordinates and the reference coordinates of the discrete nodes is calculated, and the Euclidean distance vector is decomposed into three Cartesian coordinate components, which are written into a sparse matrix corresponding to the node index to generate a geometric defect mapping tensor.

[0011] Preferably, the construction process of the real topography finite element model comprises: The initial node coordinate matrix of the standard finite element grid model is extracted, and the coordinate deviation matrix stored in the geometric defect mapping tensor is extracted; the initial node coordinate matrix and the coordinate deviation matrix are subjected to matrix linear superposition operation to generate an updated node coordinate matrix; the spatial node definition of the grid model is reset using the updated node coordinate matrix, the micro-geometric deviation is introduced under the premise of keeping the grid topology connection relationship unchanged, and a real topography finite element model is constructed.

[0012] Preferably, the process of constructing a nonlinear numerical balance equation set by a penalty function method and iteratively solving the normal interaction value between each pair of discrete contact nodes comprises: The contact major surface and the contact minor surface are defined in the real topography finite element model, and a penalty stiffness coefficient is set; the normal penetration depth of the nodes on the contact minor surface relative to the contact major surface is detected, and for the node pairs with normal penetration depth greater than zero, the penalty stiffness coefficient is added to the corresponding diagonal element position of the overall stiffness matrix; the nonlinear numerical balance equation set containing the penalty stiffness coefficient is constructed, and the Newton-Raphson algorithm is used to iteratively solve the nonlinear numerical balance equation set until the norm of the residual vector meets the preset convergence criterion; the product of the penalty stiffness coefficient and the normal penetration depth under the convergence state is output as the normal interaction value.

[0013] Preferably, the process of calculating the axial load distribution characteristic sequence based on the contact stiffness characteristic matrix comprises: constructing a global sparse diagonal matrix with the same number of nodes as the real topography finite element model; traversing all contact node pairs, filling the normal interaction value of each contact node pair as a feature element into the diagonal line position corresponding to the node index in the global sparse diagonal matrix to generate a contact stiffness feature matrix; extracting non-zero diagonal elements in the contact stiffness feature matrix, and discretizing the axial full length of the thread into a plurality of continuous axial statistical intervals; retrieving all node indexes falling into each of the axial statistical intervals, extracting the normal interaction value corresponding to each node index in the contact stiffness feature matrix, and integrating and summing the extracted normal interaction values to generate a one-dimensional load distribution feature sequence.

[0014] Preferably, the gradient convolution operation is performed on the load distribution feature sequence to calculate the local field variable peak value, specifically including: constructing a Gaussian first derivative convolution kernel, the Gaussian first derivative convolution kernel being a derivative discrete sequence of a Gaussian function, used for extracting gradient features of a signal and suppressing high-frequency noise; performing a sliding convolution operation on the load distribution feature sequence and the Gaussian first derivative convolution kernel in the time domain to output a gradient sequence reflecting the rate of change of the load distribution; searching for a modulus maximum value in the gradient sequence, locking the node position corresponding to the modulus maximum value in the time domain, and extracting the amplitude of the load distribution feature sequence at the node position as the local field variable peak value.

[0015] Preferably, the process of calculating the cumulative damage index based on the local field variable peak value and combining the damage evolution equation to output the predicted life of the target batch of high-strength bolts includes: calling an average stress correction model, substituting the local field variable peak value as an average stress parameter into the average stress correction model to calculate the fatigue life cycle number under the current load level; calculating the reciprocal of the fatigue life cycle number with respect to a single load cycle to obtain a single damage value; performing linear accumulation on the single damage values under all load conditions in the whole life cycle to generate a cumulative damage index, and determining the total cycle number when the cumulative damage index reaches a critical threshold, taking the total cycle number as the predicted life.

[0016] A fatigue life prediction system for high-strength bolt connection, comprising: a data acquisition module for acquiring a discrete point cloud coordinate sequence representing the topography of the target batch of high-strength bolts; The geometric reconstruction module is used for reconstructing the discrete point cloud coordinate sequence by a non-uniform rational B-spline interpolation method to generate a continuous geometric topology function; the preset standard thread geometric design parameter is discretized by a mesh division method to construct a standard finite element mesh model; based on the continuous geometric topology function and the standard finite element mesh model, the normal space position difference value is calculated and matrixed and packaged to generate a geometric defect mapping tensor; and based on the geometric defect mapping tensor, the node coordinate transformation is performed on the standard finite element mesh model to construct a real topography finite element model; The simulation analysis module is used for constructing a nonlinear numerical balance equation set by a penalty function method based on the real topography finite element model, iteratively solving the normal interaction value between each discrete contact node pair, and generating a contact stiffness characteristic matrix; based on the contact stiffness characteristic matrix, the load distribution characteristic sequence distributed along the axial direction is calculated, and the gradient convolution operation is performed on the load distribution characteristic sequence to calculate the local field variable peak value; The life prediction module is used for calculating the cumulative damage index based on the local field variable peak value in combination with a damage evolution equation, and outputting the predicted life of the target batch of high-strength bolts.

[0017] Compared with the prior art, the present application has the following advantages: 1、The present application constructs a three-dimensional KD tree space index to eliminate noise, and combines non-uniform rational B-spline interpolation to convert the scattered point cloud into a continuous geometric topology function; and then uses the ray equation to solve the normal deviation of the analytical surface and encapsulates it as a geometric defect mapping tensor. This substantial cooperative design solves the problem of loss of high-frequency geometric deviation (such as pitch accumulation step and flank angle disturbance) caused by smoothing processing in the prior art, ensuring the complete preservation of processing defect information in numerical simulation.

[0018] 2、The present application performs node coordinate transformation on the standard mesh through the geometric defect mapping tensor, introduces the geometric disturbance caused by the real topography under the premise of maintaining the stability of the mesh topology connection relationship; and the nonlinear numerical balance equation constructed by the penalty function method can iteratively calculate the real normal interaction between the discrete contact nodes caused by the micro topography difference. This overcomes the defect that the idealized thread model cannot calculate the local singularity of the contact stiffness matrix, so that the simulation of the local stress field and displacement field is more in line with the actual physical state.

[0019] 3、The application introduces a Gaussian first derivative convolution kernel for gradient operation after obtaining the axial load distribution sequence. This processing method can effectively suppress high-frequency noise in numerical calculation, while sharply capturing the extreme value of the change rate of the load distribution, thereby accurately positioning the stress concentration area (local field variable peak value) caused by machining errors. Combined with the damage evolution equation, the application improves the global average life prediction based on ideal geometry to local precise failure analysis based on real defects, greatly improving the robustness of the target batch high-strength bolt fatigue life prediction. BRIEF DESCRIPTION OF DRAWINGS

[0020] Figure 1 A flow chart of a high-strength bolt connection fatigue life prediction method of the application; Figure 2 A flow chart of a real topography finite element model construction of an embodiment of the application; Figure 3 A module diagram of a high-strength bolt connection fatigue life prediction system of the application. DETAILED DESCRIPTION

[0021] To make the purpose, technical scheme and advantages of the application clearer, specific embodiments will be further described below in combination with the drawings and specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the application, and not to limit the application. Other embodiments obtained by those skilled in the art based on the ideas in the present application without creative labor are also within the scope of protection of the application.

[0022] Referring to Figures 1 to 3 , the application provides a high-strength bolt connection fatigue life prediction method and system, and the specific technical scheme is as follows: Embodiment 1: A high-strength bolt connection fatigue life prediction method, referring to Figure 1 , comprising: Obtaining a discrete point cloud coordinate sequence representing the real topography of a target batch of high-strength bolts; reconstructing the discrete point cloud coordinate sequence by a non-uniform rational B-spline interpolation method to generate a continuous geometric topology function; Discretizing the pre-set standard thread geometric design parameters by a meshing method to construct a standard finite element mesh model; calculating the normal space position difference and matrix encapsulation based on the continuous geometric topology function and the standard finite element mesh model to generate a geometric defect mapping tensor; Performing node coordinate transformation on the standard finite element mesh model based on the geometric defect mapping tensor to construct a real topography finite element model; constructing a nonlinear numerical balance equation set by a penalty function method based on the real topography finite element model, and iteratively solving the normal interaction value between each discrete contact node pair to generate a contact stiffness characteristic matrix; Based on the contact stiffness feature matrix, a load distribution feature sequence distributed along the axial direction is calculated, and a gradient convolution operation is performed on the load distribution feature sequence to calculate a local field variable peak value; Based on the local field variable peak value, an accumulated damage index is calculated by combining a damage evolution equation, and a predicted service life of the target batch of high-strength bolts is output.

[0023] Further, the process of obtaining the discrete point cloud coordinate sequence representing the entity topography of the target batch of high-strength bolts includes: reading the original three-dimensional scattered point cloud data output by the optical screening device, and constructing a three-dimensional KD tree spatial index structure; traversing each data point in the original three-dimensional scattered point cloud data, searching for K neighborhood points of the data point using the three-dimensional KD tree spatial index structure, and calculating the average Euclidean distance of the data point to all K neighborhood points; calculating the mean and standard deviation of the average Euclidean distances of all data points, and removing outlier noise points whose average Euclidean distances are greater than the mean plus twice the standard deviation; mapping the remaining data points after removing the noise points from the Cartesian coordinate system to the cylindrical coordinate system, and rearranging the data points in order of increasing spiral phase angle value to generate an ordered discrete point cloud coordinate sequence.

[0024] Specifically, the original three-dimensional scattered point cloud data output by the optical screening device is read, and a three-dimensional KD tree spatial index structure is constructed. In the detection scenario of large bridge cable end anchoring bolts in this embodiment, the original three-dimensional scattered point cloud data generated by the optical screening device is a two-order tensor with a dimension of N rows and 3 columns, where N represents the total number of point clouds, usually in the order of millions, and the three columns correspond to the horizontal coordinate, vertical coordinate and depth coordinate in the three-dimensional space; the process of constructing the three-dimensional KD tree spatial index structure is as follows: calculate the numerical variance of the current point set in the three coordinate axis directions, and select the coordinate axis with the largest variance as the partition axis of the current space; extract the median value of the coordinates of all points in the partition axis direction, and use the hyperplane perpendicular to the partition axis and passing through the median value as the geometric boundary to divide the current space into two non-intersecting subspaces; recursively perform the variance calculation and space partitioning logic in each subspace until the number of data points contained in each terminal subregion reaches a preset node capacity threshold; the preset node capacity threshold is preset based on the distribution density of the point cloud sampling points, and in this embodiment, it is preset to 10, thereby establishing a hierarchical binary tree retrieval structure to provide topological support for high-speed neighbor access of large-scale point clouds.

[0025] Traverse each data point in the original three-dimensional scattered point cloud data, search K neighborhood points of the data point using the three-dimensional KD tree spatial index structure, and calculate the average Euclidean distance of the data point to all K neighborhood points. For each target point in the original data set, take its coordinate value as the query key value, perform a depth-first search in the constructed three-dimensional KD tree spatial index structure; in the search path, by comparing the target point coordinates with the distance of each level partition plane, exclude the spatial branches that do not contain potential neighbor points, until the K points closest to the spatial position of the target point are locked, K in the K neighborhood points is the preset number of neighbor search, and in this embodiment, K is preset to 20 according to the scanning point density per unit area of the bolt surface; calculate the spatial Euclidean distance of the target point to each neighborhood point, and the calculation logic is: calculate the coordinate difference of two points in three coordinate axis directions, perform square operation on the three difference values respectively and sum, and then perform arithmetic square root operation on the sum result; accumulate the obtained 20 spatial Euclidean distances, and divide the accumulated sum by 20 to obtain the average Euclidean distance of the target point, which represents the local discrete degree of the point on the bolt surface.

[0026] Calculate the mean and standard deviation of the average Euclidean distances of all data points, and remove outlier noise points whose average Euclidean distance is greater than the mean plus twice the standard deviation. By performing an arithmetic average operation on the average Euclidean distances of all points in the data set, the mean value reflecting the overall density of the point cloud is obtained; calculate the square of the difference between the average Euclidean distance of each point and the mean value, sum all square values and divide by the total number of points to perform an arithmetic square root operation to obtain the standard deviation; construct a noise judgment rule based on statistical boundaries: if the average Euclidean distance of a data point is greater than the sum of the mean value and twice the standard deviation, it is determined that the point is an outlier noise point due to strong reflection of the metal surface or shielding of dust, and the point is deleted from the tensor data; the threshold of twice the standard deviation is set according to the statistical characteristics of normal distribution, which can retain the fine geometric features of the thread root while removing isolated measurement noise points deviating from the bolt entity.

[0027] The remaining data points after the noise points are removed are mapped from the Cartesian coordinate system to the cylindrical coordinate system, and the data points are rearranged in order of increasing spiral phase angle value to generate an ordered discrete point cloud coordinate sequence. Before performing the mapping, first, the principal component analysis algorithm is used to determine the central principal axis vector of all point sets, and the vector is defined as the longitudinal height axis of the cylindrical coordinate system; for each remaining data point, the vertical distance of the projection point of the data point in the plane perpendicular to the longitudinal height axis to the principal axis is calculated, which is defined as the radial distance; the angle of the projection point in the plane relative to the preset reference axis is calculated, which is obtained by calculating the inverse tangent value of the ratio of the longitudinal coordinate to the transverse coordinate, and is defined as the spiral phase angle, which is a physical quantity describing the circumferential position of the point on the spiral thread path; considering the continuous growth characteristics of the bolt thread, the periodic jump of the phase angle at the closed loop position of the circumference is detected, and the change amount of the longitudinal height is combined to calculate the cumulative phase angle value of each point on the full spiral path; the cumulative phase angle values are sorted in ascending order from small to large, so that the originally randomly distributed data points in three-dimensional space are converted into an ordered discrete point cloud coordinate sequence with M rows and 3 columns (M is the number of points after denoising) and strictly following the physical rotation direction logic of the thread.

[0028] Through the synergistic effect of the three-dimensional KD tree and the statistical distance determination algorithm, automatic filtering and spatial ordering of the bridge anchoring bolt surface measurement noise are realized. This processing procedure converts the unstructured geometric measurement tensor into an ordered sequence that conforms to the thread topology logic, solving the technical pain point that discrete scanning data is difficult to directly align with finite element nodes due to disorder, and providing a high-fidelity and structured data carrier for subsequent establishment of a continuous geometric topology function taking into account micro defects.

[0029] Further, the generation process of the continuous geometric topology function includes: calculating the parameter node vector of each data point in the discrete point cloud coordinate sequence using the cumulative chord length parameterization method; constructing a linear approximation equation set based on the least squares method, substituting the discrete point cloud coordinate sequence as the target observation value into the linear approximation equation set, and solving the control vertex coordinates that minimize the equation set residual sum of squares; using the control vertex coordinates, the parameter node vector, and the preset basis function order, constructing a parameterized continuous geometric topology function, which is an analytical expression describing the geometric trajectory of the thread surface.

[0030] Specifically, the parameter node vector of each data point in the discrete point cloud coordinate sequence is calculated using the cumulative chord length parameterization method. In the micro-topography reconstruction process of the bridge anchoring bolt in this embodiment, since the bolt thread presents a spiral ascending spatial form, and the sampling density of discrete scanning points is not completely uniform at local features, direct use of equal interval parameterization will cause severe oscillation or geometric distortion in surface fitting.

[0031] The cumulative chord length parameterization method refers to a calculation logic for assigning a corresponding normalized parameter value to each data point in a discrete point cloud coordinate sequence based on the path arc length weight between each measurement point in a physical space. In the embodiment, the specific steps performed are as follows: calculating the spatial straight line distance between each pair of adjacent data points in an ordered sequence of discrete point cloud coordinates, i.e., calculating the coordinate difference of two adjacent points in the three coordinate axis directions of the Cartesian coordinate system, performing square operation on the three difference values respectively and summing them up, and then performing arithmetic square root operation on the sum to obtain the chord length value of each segment; calculating the cumulative sum of all adjacent segment chord length values from the starting point of the ordered sequence to the currently calculated data point; and simultaneously calculating the total cumulative length of all adjacent segment chord length values in the entire sequence.

[0032] Dividing the cumulative sum corresponding to each data point by the total cumulative length obtains a numerical sequence distributed in the closed interval of zero to one and showing a monotonically increasing trend, i.e., a parameter node vector, which reflects the relative physical positions of the points on the surface of the bridge bolt on the real spiral trajectory; through this method, even if the point cloud is densely distributed at the thread tip or root, the parameter interval assigned to it will be correspondingly increased or decreased, thereby ensuring that the generated parameter distribution is strictly proportional to the span of the physical appearance of the bolt, and providing a precise parameter benchmark for subsequent construction of a lossless analytical surface.

[0033] A linear approximation equation set based on the least squares method is constructed, and the discrete point cloud coordinate sequence is substituted into the linear approximation equation set as the target observation value to solve the control vertex coordinates that minimize the residual sum of squares of the equation set. The control vertex coordinates refer to a set of feature points used to define the shape of the non-uniform rational B-spline surface, the number of which is pre-set to a fixed value less than the total number of discrete point clouds, and in the embodiment, the number of control vertices is pre-set to ten percent of the total number of point clouds to achieve feature compression of massive data. The construction logic of the linear approximation equation set is as follows: taking the pre-set basis function order and the aforementioned parameter node vector as known independent variables, calculating the non-uniform rational B-spline basis function matrix corresponding to each parameter node as the coefficient matrix of the linear equation set; taking the input discrete point cloud coordinate sequence as the target observation vector; and the solving process follows the least squares criterion, i.e., constructing a normal equation set through matrix transposition and matrix multiplication operation, and calculating a set of optimal control vertex coordinate matrices that minimize the distance square sum between the predicted space coordinates of the fitted surface at each parameter node and the actual discrete measurement coordinates.

[0034] The parameterized continuous geometric topology function is constructed by using the control vertex coordinates, the parameter node vector and a preset base function order. The preset base function order is a numerical value set based on the continuity requirement of the thread surface, and is preset to cubic (third order) in this embodiment, so as to ensure that the generated surface has a second-order continuous derivative. The specific construction logic is that the calculated control vertex coordinates are subjected to linear weighted summation with the corresponding third-order non-uniform rational B-spline base function. The obtained continuous geometric topology function is an analytical mapping relationship with the parameter node vector as an input variable and the three-dimensional space coordinates as an output variable. The function not only can represent the real pitch fluctuation and surface texture of the bridge bolt, but also can smoothly connect the transition area of the tooth base and the tooth side due to its mathematical continuity, so as to convert the discrete measurement point traces into an analytical surface with physical meaning and differentiability.

[0035] The accurate alignment of the physical path and the parameter domain is realized by the accumulated chord length parameterization method, and the geometric distortion caused by uneven sampling density is eliminated. In combination with the least square linear approximation, the embodiment realizes high-fidelity reconstruction of the micro-machining defects of the bridge bolt while effectively compressing the massive detection data. The generated continuous geometric topology function has analytical continuity, solves the problem that the discrete point cloud is difficult to calculate the normal vector due to the lack of differential characteristics, and provides a reliable geometric source for subsequent high-precision solution of the contact stiffness matrix.

[0036] Further, with reference to Figure 2 , the generation process of the geometric defect mapping tensor includes: reading the standard thread geometric design parameters to generate an ideal spiral surface, and discretizing the ideal spiral surface in space by using a hexahedral element to construct a standard finite element grid model including a node coordinate matrix and an element connection relationship; traversing the discrete nodes located on the contact surface in the standard finite element grid model, and constructing a ray equation along the unit normal vector of each discrete node; simultaneously solving the ray equation and the continuous geometric topology function, and solving the spatial intersection point coordinates of the analytical surface defined by the ray equation and the continuous geometric topology function by using a numerical root solving algorithm; calculating the Euclidean distance vector between the spatial intersection point coordinates and the reference coordinates of the discrete nodes, and decomposing the Euclidean distance vector into three Cartesian coordinate components, writing the Cartesian coordinate components into a sparse matrix corresponding to the node index, and generating a geometric defect mapping tensor.

[0037] Specifically, the standard thread geometry design parameters are read to generate an ideal helical curved surface, and hexahedral elements are used to discretize the ideal helical curved surface in space to construct a standard finite element grid model including a node coordinate matrix and an element connection relationship. In the numerical simulation scene of the bridge anchoring bolt in the embodiment, the generation logic of the ideal helical curved surface is as follows: first, according to the preset nominal diameter (such as 100 mm), pitch (such as 6 mm) and thread angle (such as 60 degrees), the parameter equation about the radial distance, phase angle and axial height is established in the cylindrical coordinate system. The parameter equation describes the synchronous rotation and translation motion of a cross section conforming to the 60-degree equilateral triangle profile along the central axis, thereby forming an analytically continuous helical envelope in three-dimensional space.

[0038] The hexahedral element refers to a three-dimensional entity isoparametric element with eight corner points, six planes and twelve edges. The linear hexahedral element with reduced integration technology (i.e. C3D8R element) is used in the embodiment, which sets an integration point at the center of the element to calculate stress and strain, so as to avoid the shear self-locking phenomenon in the thread tooth stress analysis.

[0039] The specific logic of discretizing the ideal helical curved surface in space and constructing the standard finite element grid is as follows: the sweep division algorithm is used to divide the radial thickness direction of the thread entity into 5-10 layers of equal interval segments, the circumferential phase angle is divided into 120-200 angular segments according to 360 degrees per circle, and the axial height is divided into several slices according to the pitch; the spatial nodes are defined at the intersection points in the radial, circumferential and axial directions to generate a node coordinate matrix composed of node numbers and their horizontal coordinates, vertical coordinates and depth coordinates in the Cartesian coordinate system; at the same time, the eight node numbers associated with each element are defined in counterclockwise order to form an element connection relationship table. The dimension of the node coordinate matrix is the total number of grid nodes in rows and 3 in columns; the dimension of the element connection relationship table is the total number of elements in rows and 8 in columns. The model is used as a pure geometric reference model for subsequent processing defects.

[0040] The discrete nodes on the contact surface in the standard finite element mesh model are traversed, and a ray equation is constructed along the unit normal vector of each discrete node. For the nodes on the force side of the thread tooth in the standard finite element mesh model, node normal extraction is performed: for each target node, a unit of element connection relationship table is retrieved to determine all surface element patches that are co-vertices with the target node; the unit normal vector of each triangle or quadrilateral patch is calculated; the unit normal vectors of all adjacent patches are subjected to an arithmetic average operation to obtain the average unit normal vector of the node; the ideal position coordinates of the node are taken as the starting vector, the average unit normal vector is taken as the direction vector, and a to-be-determined length scalar parameter is introduced to construct a ray equation describing a straight line trajectory in space. The equation represents a mathematical path that starts from the ideal geometric point and searches for the intersection point of the real physical surface perpendicular to the theoretical surface.

[0041] The ray equation and the continuous geometric topology function are combined, and the spatial intersection coordinates of the analytical surface defined by the ray equation and the continuous geometric topology function are solved by a numerical root-finding algorithm; when there are multiple spatial intersection points, the Euclidean distances between each intersection point and the original discrete node are calculated, and the spatial intersection coordinates of the effective spatial intersection point with the smallest Euclidean distance are selected to ensure the uniqueness of the mapping. The numerical root-finding algorithm is a nonlinear equation approximation solving technique based on the Newton-Raphson iteration logic. In the intersection calculation of the real topography surface of the bridge bolt, the specific execution logic is as follows: the three-dimensional space coordinate expression in the ray equation is substituted into the continuous geometric topology function based on the non-uniform rational B-spline constructed in the foregoing step, a nonlinear residual function with the ray length scalar parameter as the only unknown quantity is constructed, and the function represents the spatial distance between a point on the ray and the analytical surface. In order to ensure the uniqueness of the solution under complex micro-topography, the initial value of the length scalar parameter is set to zero. An iterative update loop is executed: the residual value under the current scalar parameter is calculated; the derivative of the residual function is calculated to calculate the rate of change of the residual with respect to the scalar parameter; the residual value is divided by the rate of change to obtain the correction increment of the scalar parameter; the correction increment is subtracted from the current scalar parameter to complete one parameter update; the loop logic is repeated until the absolute value of the scalar parameter change between two consecutive updates is less than the preset convergence threshold, and the preset convergence threshold is set to 0.001 microns based on the bridge bolt machining precision requirement. The final converged length scalar parameter is substituted back into the ray equation to calculate the accurate intersection coordinates of the ray and the analytical surface defined by the continuous geometric topology function in three-dimensional space. By adopting the zero initial value strategy, the algorithm is constrained to search within the ideal surface neighborhood, thereby ensuring that the calculated intersection coordinates are the physically effective solution with the smallest Euclidean distance, and the mathematical multiple solution ambiguity possibly caused by surface distortion is excluded.

[0042] The Euclidean distance vector between the spatial intersection point coordinates and the reference coordinates of the discrete nodes is calculated, and the Euclidean distance vector is decomposed into three Cartesian coordinate components, and written into a sparse matrix corresponding to the node index, to generate a geometric defect mapping tensor. The geometric defect mapping tensor is a structured data tensor for storing the micro-topography deviation values of the bridge bolts. In order to ensure that the node definition of the finite element model can be directly operated algebraically, the tensor is constructed as a real global sparse matrix with the number of rows equal to the total number of standard finite element grid nodes and the number of columns equal to three in specific implementation. The calculation logic is as follows: the three-dimensional values of the spatial intersection point coordinates obtained by solving are subtracted from the corresponding reference coordinates of the discrete nodes respectively to obtain a three-dimensional Euclidean distance vector. The vector represents the normal displacement of the real surface relative to the ideal design surface caused by the processing error in the physical sense. Then, the displacement vector is decomposed into the horizontal axis component, the vertical axis component and the vertical axis component in the Cartesian coordinate system, and is sequentially filled into the row position and the corresponding column position determined by the current node global index number in the sparse matrix. For the nodes of the remaining non-contact surface in the grid model, the matrix row element values corresponding to the nodes are all preset to zero. The generated geometric defect mapping tensor can directly participate in the operation as a coordinate increment item by establishing the same data dimension as the initial node coordinate matrix, so as to realize the digital integration of the processing quality data in the finite element model.

[0043] By constructing a numerical root-finding mechanism based on the normal ray tracing, the accurate spatial alignment between the ideal finite element grid and the discrete measurement surface is realized. The processing method encapsulates the micro-processing deviation as a standardized geometric defect mapping tensor while keeping the grid topological connection relationship completely consistent. This solves the problem that the existing discrete point cloud cannot directly drive the deformation of the continuous grid, and ensures that the micro-geometric defects of the bridge bolt can be transmitted to the numerical balance equation in the form of node coordinate correction without loss, providing accurate geometric boundary constraints for subsequent solving of the non-uniform distribution of local contact stiffness.

[0044] Further, with reference to Figure 2 , the construction process of the real topography finite element model includes: extracting the initial node coordinate matrix of the standard finite element grid model, and extracting the coordinate deviation matrix stored in the geometric defect mapping tensor; performing matrix linear superposition operation on the initial node coordinate matrix and the coordinate deviation matrix to generate an updated node coordinate matrix; resetting the spatial node definition of the grid model by using the updated node coordinate matrix, introducing the micro-geometric deviation under the premise of keeping the grid topological connection relationship unchanged, and constructing the real topography finite element model.

[0045] Specifically, an initial node coordinate matrix of the standard finite element mesh model is extracted, and a coordinate deviation matrix stored in the geometric defect mapping tensor is extracted. In the simulation model reconstruction of the bridge anchoring bolt in the embodiment, the initial node coordinate matrix is a second-order tensor derived from the aforementioned standard finite element mesh model, and the dimension is the total number of mesh nodes multiplied by 3. Each row corresponds to the index number of a node, and the three columns store the horizontal coordinate, the vertical coordinate and the axial height of the node in the ideal spiral geometric state. The coordinate deviation matrix is an effective numerical set extracted from the sparse geometric defect mapping tensor, and the dimension is completely consistent with that of the initial node coordinate matrix. In the coordinate deviation matrix, only the node rows belonging to the thread contact surface store non-zero offset values, and these values represent the normal deviation of the actual profile from the ideal profile obtained by optical measurement and numerical root solving. The node rows of the non-contact area (such as the non-mating area of the bolt shaft) are preset to zero to ensure the local pertinence of the model reconstruction.

[0046] The initial node coordinate matrix and the coordinate deviation matrix are subjected to linear superposition operation to generate an updated node coordinate matrix. The specific logic for generating the updated node coordinate matrix is as follows: first, a data dimension verification process is performed to ensure that the initial node coordinate matrix and the aforementioned generated geometric defect mapping matrix are strictly consistent in data dimension, i.e., both are real matrices with the number of rows equal to the total number of mesh nodes and the number of columns equal to three; then, a linear superposition operation based on matrix algebra is performed: traverse each row of the matrix, extract the horizontal axis component, the vertical axis component and the vertical axis component stored in the current row of the geometric defect mapping matrix, and perform arithmetic summation with the ideal horizontal coordinate, the ideal vertical coordinate and the ideal axial height corresponding to the same row in the initial node coordinate matrix; by performing the above element-by-element addition operation on the full matrix, the discrete geometric displacement deviation is directly compensated to each ideal node, thereby synthesizing an updated node coordinate matrix with the dimension of the total number of nodes multiplied by 3, but the numerical value has been fused with the bridge bolt entity processing deviation (such as the non-uniform increment of the pitch); finally, the spatial node definition of the finite element mesh model is reset using the updated node coordinate matrix, thereby successfully constructing a real profile finite element model while keeping the original mesh topology connection relationship completely unchanged.

[0047] The space node definition of the grid model is reset by using the updated node coordinate matrix, micro-geometric deviations are introduced while the grid topology connection relationship is maintained, and a real topography finite element model is constructed. The specific logic of constructing the real topography finite element model is as follows: the preset unit connection relationship table in the standard finite element grid model is maintained, which records the mapping logic of the number of each hexahedral unit and its corresponding eight nodes; the replacement and reset of the node space position is performed, that is, the coordinate values pointed to by each node number in the simulation calculation domain are switched from the ideal values in the initial node coordinate matrix to the real values in the updated node coordinate matrix; this process realizes the controlled deformation of the grid skeleton at the mathematical bottom layer, which does not change the total number of units and the adjacency relationship between nodes of the grid, but changes the space position of the nodes, so that the originally regular hexahedral unit generates micro-geometric distortion conforming to the physical entity defect in the tooth side surface area; finally, the generated real topography finite element model has the ability to describe the micro-topography disturbance and the locally singular contact stiffness caused by it at the calculation level.

[0048] As a preferred embodiment, the construction process of the real topography finite element model includes: calculating the Jacobian determinant value of each hexahedral unit after node coordinate transformation is performed; retrieving an abnormal unit set whose Jacobian determinant value is less than a preset safety threshold, and determining the disturbed node index involved in the abnormal unit set; calling the geometric defect mapping tensor, performing Laplacian smoothing filtering on the coordinate deviation value corresponding to the disturbed node index involved, and generating a corrected mapping tensor; re-executing the node coordinate transformation based on the corrected mapping tensor until the Jacobian determinant values of all model units are greater than the preset safety threshold, and the model construction is completed.

[0049] Specifically, the Jacobian determinant value of each hexahedral unit after the node coordinate transformation is performed is calculated. In a large-scale hexahedral grid for a bridge anchoring bolt, the Jacobian determinant value represents the volume linear transformation ratio of the unit from the standard parametric space to the physical space. The calculation logic is as follows: the eight vertex coordinates of each unit are extracted, the partial derivative matrix of the unit geometric function with respect to the local natural coordinate system is calculated, and the determinant operation is performed on the matrix. If the value is close to zero or negative, it is determined that the unit has collapsed in physical sense due to the introduction of excessive geometric disturbance.

[0050] The abnormal unit set is retrieved and the disturbed nodes are locked. The preset safety threshold is an empirical value based on the convergence robustness of the numerical solver, which is set to 0.15 in this embodiment. The system traverses all the units of the whole model, and locks the unit indexes whose Jacobian numerical value is lower than 0.15. Then, by looking up the unit connection relationship table, the node numbers contained in these distorted units are determined and recorded as the disturbed nodes. These nodes are usually located at the micro positions where the surface processing burr of the bridge bolt is most sharp or the point cloud scanning noise is left.

[0051] The Laplace smoothing filter based on the geometry defect mapping tensor is performed. Specifically, for each disturbed node, all its physically adjacent nodes are searched by using the grid topology connection relationship; the coordinate deviation components stored in the geometry defect mapping tensor of these neighbor nodes are extracted; the arithmetic mean of the deviations of all neighbor nodes is calculated, and the mean value is used to replace the original deviation value of the current disturbed node, thereby generating a corrected mapping tensor. This process reduces the local gradient that is sufficient to destroy the stability of the grid topology while preserving the main features of the processing deviation, achieving a balance between physical reality and grid quality.

[0052] By introducing the unit Jacobian determinant constraint, the decoupling between the processing defect injection and the numerical calculation stability is realized, ensuring that the model still has high solving convergence after taking into account the complex processing deviation.

[0053] By performing linear superposition of the node coordinate matrix and the coordinate deviation matrix, the lossless injection of the processing deviation information is realized without changing the finite element grid topology structure. This processing method avoids the waste of computing resources and data alignment errors caused by re-dividing the grid, ensuring that the simulation model can accurately bear the micro-geometric defects on the surface of the bridge anchor bolt. The generated real topography finite element model can truly reflect the tooth side topography disturbance caused by uneven processing, providing a high-fidelity geometric reference for accurately identifying the local field variable peak value and evaluating the fatigue danger zone.

[0054] Further, the process of constructing a nonlinear numerical balance equation set by a penalty function method and iteratively solving the normal interaction value between each pair of discrete contact nodes includes: defining a contact major surface and a contact minor surface in the real topography finite element model and setting a penalty stiffness coefficient; detecting the normal penetration depth of the nodes on the contact minor surface relative to the contact major surface, and for the node pairs with a normal penetration depth greater than zero, adding the penalty stiffness coefficient to the corresponding diagonal element position of the overall stiffness matrix; constructing the nonlinear numerical balance equation set containing the penalty stiffness coefficient, and iteratively solving the nonlinear numerical balance equation set by using the Newton-Raphson algorithm until the norm of the residual vector meets the preset convergence criterion; and outputting the product of the penalty stiffness coefficient and the normal penetration depth in the convergence state as the normal interaction value.

[0055] Specifically, the contact master surface and the contact slave surface are defined in the real topography finite element model, and a penalty stiffness coefficient is set. In the simulation analysis of the bridge anchoring bolt system, the contact master surface is defined on the nut inner thread surface with higher structural stiffness, and the geometric topology is kept in an analytical ideal state; the contact slave surface is defined on the bolt outer thread surface with micro machining defects, and the node space definition has been updated according to the geometric defect mapping tensor. The penalty stiffness coefficient refers to a virtual elastic coefficient used in the penalty function contact algorithm to generate numerical resistance at the contact interface to limit geometric penetration. The specific logic of setting the penalty stiffness coefficient is: extracting the preset elastic modulus of the bridge bolt material (such as the elastic modulus of alloy steel is 210000 MPa) and the characteristic element size of the grid model (such as 0.2 mm); calculate the ratio of the elastic modulus and the characteristic element size as the reference stiffness. The preset penalty stiffness coefficient is determined by performing multiple trials on the validation set, and the setting principle is to make the maximum penetration depth less than one thousandth of the characteristic element size on the premise of ensuring the convergence speed of the numerical solution. In this embodiment, the penalty stiffness coefficient is preset to 100 times the reference stiffness, and the specific numerical order is set to 10 13 Newton per cubic meter. This setting balances between physical reality and numerical stability, ensuring that the contact stress caused by the micro concave-convex of the thread side can be stably captured.

[0056] The normal penetration depth of the nodes on the contact slave surface relative to the contact master surface is detected, and for the node pairs with normal penetration depth greater than zero, the penalty stiffness coefficient is added to the corresponding diagonal element position of the overall stiffness matrix. The normal penetration depth refers to the normal displacement of the contact slave surface node into the contact master surface geometric boundary under the action of load. The detection logic is: in each increment step, calculate the shortest distance from the slave surface node to the discrete topology surface of the master surface, if the distance direction is opposite to the normal direction of the master surface, it is determined that penetration occurs; for the node pairs that occur penetration, the preset penalty stiffness coefficient is directly added to the diagonal element corresponding to the node degree of freedom in the overall stiffness matrix as an additional stiffness term; the overall stiffness matrix is a symmetric sparse matrix with a dimension of the total number of model degrees of freedom multiplied by the total number of model degrees of freedom. This processing process is equivalent to introducing a fictitious spring with extremely high stiffness at the micro topography position where interference occurs, to prevent geometric overlap between physical entities.

[0057] As a preferred embodiment, the iterative solution of the normal interaction value between each pair of discrete contact nodes comprises: extracting the two principal curvature values of the continuous geometric topology function at the physical position corresponding to each pair of discrete contact nodes; calculating the geometric mean of the two principal curvature values to generate a local topographic complexity factor; constructing a penalty stiffness adjustment weight based on the local topographic complexity factor, and performing hierarchical correction on the preset penalty stiffness coefficient using the penalty stiffness adjustment weight to generate an adaptive penalty stiffness distribution sequence; and constructing a nonlinear numerical balance equation set using the adaptive penalty stiffness distribution sequence.

[0058] Specifically, the principal curvature values of the continuous geometric topology function at the position corresponding to each pair of discrete contact nodes are extracted. The calculation logic is: calling the analytical expression system constructed by the non-uniform rational B-spline, performing second-order partial derivative operation on the surface coordinates with respect to the parameter space, and constructing the first and second fundamental form matrices of the surface geometry; solving the matrix characteristic equation to obtain the two principal curvatures representing the maximum and minimum bending degree at the position. This value accurately quantifies the "steepness" of the bridge bolt surface caused by processing at the micro level.

[0059] The local topographic complexity factor is generated and the penalty stiffness adjustment weight is constructed. The calculation logic is: performing arithmetic multiplication on the absolute values of the two principal curvature values and taking the square root to obtain the geometric mean scalar of the Gaussian curvature, which is defined as the local topographic complexity factor. Based on this factor, the adjustment weight is set: when the complexity factor is high (representing micro convex peaks), the adjustment weight is preset to a proportionality coefficient between 1.2 and 1.5; when the complexity factor is low, the adjustment weight is preset to 1.0. The setting is based on the fact that the effective contact stiffness at the micro convex position is significantly higher than the average value due to stress shrinkage effect.

[0060] The adaptive penalty stiffness distribution sequence is generated and the balance equation set is constructed. The execution logic is: the preset penalty stiffness coefficient (such as 100 times the material elastic modulus) is compounded with the adjustment weight corresponding to each contact node one by one. The generated distribution sequence is a numerical vector that fluctuates with the spatial topography. When constructing the nonlinear numerical balance equation set, the uniform penalty term in the diagonal element position of the global stiffness matrix is replaced by the adaptive distribution sequence corresponding to its spatial index, so that the solver can perceive the non-uniformity of the contact stiffness caused by the micro topography.

[0061] By introducing the adaptive penalty stiffness adjustment mechanism based on principal curvature, the physical coupling of contact force solution and micro topographic features is realized, the underestimation error of traditional algorithms for micro convex point pressure is greatly corrected, and the physical fidelity of stress field solution is improved.

[0062] The nonlinear numerical balance equation set containing the penalty stiffness coefficient is constructed, and the Newton-Raphson algorithm is used to iteratively solve the nonlinear numerical balance equation set until the norm of the residual vector meets the preset convergence criterion. The nonlinear numerical balance equation set refers to a set of algebraic equations describing the internal force, external force and contact constraint force in the bolted connection system to achieve static equilibrium state. The specific process of construction is as follows: first, according to the node displacement state of the real topography finite element model, the node internal force vector is calculated, which is obtained by integrating the element stress field in the element volume and represents the elastic resistance of the material itself; second, the externally applied bridge pre-tightening load is extracted to construct the external load vector; third, the contact force vector is constructed, and for each pair of contact nodes with a normal penetration depth greater than zero, the contact force component is defined as the product of the preset penalty stiffness coefficient and the current normal penetration depth; finally, the external load vector is subtracted from the node internal force vector, and then the contact force vector generated by the penalty function is subtracted, and the difference is set to zero, thereby forming a nonlinear numerical balance equation set with node displacement as the independent variable. In this equation set, the penalty stiffness coefficient directly participates in the construction of the tangent stiffness matrix of the equation set as a contact stiffness term, so that the equation set can perceive the contact state switching caused by micro defects.

[0063] The execution process of the Newton-Raphson algorithm is as follows: in the current incremental step, the overall tangential stiffness matrix containing the penalty stiffness contribution is calculated using the current node position; the current residual vector is calculated, which is the difference between the external load vector and the sum of the current internal force and contact force vectors; the linearized equation set is solved to obtain the modified increment of node displacement; the spatial coordinates of the nodes are updated using the modified increment, and the penetration detection and stiffness superposition logic are re-executed. Repeat the above iteration process until the norm of the residual vector (i.e. the square root of the sum of the squares of the unbalanced components of force) is less than the preset convergence criterion. The preset convergence criterion is set based on the calculation accuracy requirement that the residual norm is less than 0.5% of the average contact force in the current load step, ensuring that a high-precision static equilibrium solution is obtained under complex geometric constraints containing micro defects.

[0064] The product of the penalty stiffness coefficient and the normal penetration depth in the output convergence state is taken as a normal interaction value. The normal interaction value refers to the normal contact force transmitted between the contact node pairs in the equilibrium state. The calculation logic is as follows: after the Newton-Raphson algorithm iteration ends and the convergence condition is reached, the final residual normal penetration depth value of each contact node is extracted, and is subjected to an arithmetic multiplication operation with a preset penalty stiffness coefficient; the obtained value can accurately reflect the real load transmission strength of each discrete contact position of the bridge bolt under the real surface topography containing micro defects. Due to the introduction of the real machining deviation in the model, the obtained normal interaction value will present significant non-uniformity in spatial distribution, and can represent the local stress concentration phenomenon caused by pitch fluctuation or flank angle disturbance.

[0065] By introducing the penalty function method and the Newton-Raphson iteration mechanism, the nonlinear contact balance problem considering micro defects is accurately solved. The method can effectively capture the local singularity of contact stiffness caused by the machining deviation of the bridge bolt, and convert the abstract geometric defects into specific normal interaction force distribution. This overcomes the defect that the contact force distribution is too smooth under the ideal geometric model, and provides real physical feedback for subsequent accurate extraction of local field variable peak value and evaluation of fatigue danger zone.

[0066] Further, the process of calculating the load distribution feature sequence distributed along the axial direction based on the contact stiffness feature matrix comprises: constructing a global sparse diagonal matrix with a dimension equal to the total number of nodes of the real topography finite element model; traversing all contact node pairs, taking the normal interaction value of each contact node pair as a feature element, filling it into the diagonal line position corresponding to the node index in the global sparse diagonal matrix, and generating a contact stiffness feature matrix; extracting the non-zero diagonal elements in the contact stiffness feature matrix, and discretizing the axial full length of the thread into a plurality of continuous axial statistical intervals; retrieving all node indexes falling into each of the axial statistical intervals, extracting the normal interaction values corresponding to the node indexes in the contact stiffness feature matrix, and integrating and summing the extracted normal interaction values to generate a one-dimensional load distribution feature sequence.

[0067] Specifically, a global sparse diagonal matrix with dimensions equal to the total number of nodes in the actual morphology finite element model is constructed. This global sparse diagonal matrix is ​​a square matrix with both the number of rows and columns equal to the total number of mesh nodes N. This matrix has non-zero storage space only on the main diagonal. The specific construction logic is as follows: First, the updated node coordinate matrix stored in the actual morphology finite element model is retrieved, and its row number is obtained as the total number of nodes N. Then, memory space is allocated based on the total number of nodes N to construct a second-order square matrix data structure, and all elements in this square matrix are initialized to zero. To optimize computational efficiency, this matrix adopts a sparse storage format, that is, physical storage addresses are allocated only to diagonal elements with equal row and column indices, thereby establishing a numerical carrier that can map discrete mechanical responses one-to-one with the global mesh topology.

[0068] Iterate through all contact node pairs, using the normal interaction values ​​of each contact node pair as feature elements, and fill them into the diagonal positions corresponding to the node indices in the global sparse diagonal matrix to generate a contact stiffness feature matrix. The specific assembly logic for generating the contact stiffness feature matrix is ​​as follows: First, call the contact force data table obtained from solving the aforementioned nonlinear numerical equilibrium equations, and filter out the list of all face node indices with normal penetration depths greater than zero; for each node index number k in the list, extract its corresponding normal interaction value. The normal interaction value is obtained by calculating the product of a preset penalty stiffness coefficient and the normal penetration depth of the node in the convergent state. Subsequently, the matrix element filling assignment operation is performed: using the index number k as the addressing key of the matrix, the intersection of the kth row and the kth column of the global sparse diagonal matrix is ​​located, and the extracted normal interaction value is written into the cell. By performing the above "point-to-point" numerical filling on all nodes in the contact state, a contact stiffness feature matrix is ​​generated. In this embodiment, the contact stiffness feature matrix is ​​defined as a data carrier used to characterize the degree of response of penalty stiffness at micro-defects. Although the physical dimension of the numerical sequence stored on its diagonal is force, its numerical distribution directly reflects the non-uniform state of contact stiffness caused by morphological differences. This matrix intuitively quantifies the difference in local contact bearing capacity of bridge bolts caused by processing morphological defects through the spatial fluctuation of numerical values, and transforms the discrete response originally implicit in the contact algorithm into a global feature vector that can perform matrix operations.

[0069] The non-zero diagonal elements in the contact stiffness feature matrix are extracted, and the axial length of the thread is discretized into multiple continuous axial statistical intervals. When processing the macroscopic load distribution of bridge bolts, it is necessary to convert high-dimensional surface contact features into one-dimensional signals distributed along the bolt axis. This conversion process relies on the spatial discretization of the axial length. The specific execution logic of discretization is as follows: obtain the maximum and minimum values ​​of the axial coordinates in the updated node coordinate matrix, calculate their difference to determine the effective axial length of the bolt thread segment (e.g., 500 mm). A preset interval step size characterizing the spatial sampling resolution is used. In this embodiment, the preset step size is 1.0 mm, based on one-sixth of the pitch. This value is set to ensure that each sampling interval contains complete tooth profile features to suppress local mesh fluctuations. The axial length is divided into multiple continuous axial statistical intervals according to the interval step size. Each axial statistical interval is defined as a thin spatial layer with a defined start and end range of axial coordinates.

[0070] Retrieve all node indices falling within each of the axial statistical intervals, extract the normal interaction values ​​corresponding to each node index in the contact stiffness feature matrix. Considering that fatigue failure usually originates at the local location most severely affected by micro-defects rather than the overall cross-section under stress, feature filtering is performed on the normal interaction values ​​of all nodes in the interval, extracting the top 5% of nodes with the largest values ​​and calculating their arithmetic mean, which is used as the characteristic contact strength at that axial location, thereby generating a one-dimensional load distribution feature sequence. The core of the transformation process lies in spatial alignment and local extreme value feature extraction: Each axial statistical interval is traversed, and all nodes whose current physical spatial position falls within that interval are selected based on the axial coordinate values ​​(Z coordinates) in the node coordinate matrix; the normal interaction values ​​corresponding to the selected node indices are extracted from the diagonal of the contact stiffness feature matrix; considering that fatigue failure usually originates in the area of ​​highest local stress concentration rather than the overall cross-section under stress, feature extraction is performed on the extracted values: the normal interaction values ​​of all nodes within the interval are sorted by value, and the top 5% of nodes with the largest values ​​are selected. The arithmetic mean of the normal interaction values ​​of these nodes is calculated and defined as the characteristic contact strength at that axial position. By repeating the above selection and feature extraction process for all axial statistical intervals, a numerical sequence composed of the characteristic contact strength of each interval is finally generated, i.e., a one-dimensional load distribution feature sequence. This sequence is a vector with one row and one column, and its fluctuation trend accurately reflects the stress bearing change trajectory of bridge bolts at local danger points along the axial direction, even with microscopic defects.

[0071] The process of generating a one-dimensional load distribution feature sequence also includes performing numerical normalization. This numerical normalization refers to the computational logic of mapping the characteristic contact intensity of each statistical interval to a unified dimension space. Specifically, the process involves calculating the arithmetic mean of the characteristic contact intensity of all axial statistical intervals, and then dividing the characteristic contact intensity of each interval by this arithmetic mean to eliminate the influence of numerical sensitivity caused by differences in the magnitude of external preload. The normalized interval values ​​are then vectorized and arranged in order of axial position from nearest to farthest, ultimately generating a one-dimensional load distribution feature sequence. This sequence, in the form of a dimensionless fluctuating signal, accurately characterizes the relative distribution characteristics of local peak values ​​of axial loads caused by microscopic defects in bridge bolts.

[0072] By constructing a global sparse diagonal matrix and performing axial interval-based integral summation, dimensionality reduction extraction of complex local contact stress fields and macroscopic load distribution characteristics is achieved. This processing method effectively eliminates high-frequency numerical interference caused by mesh discretization. Furthermore, by axially reorganizing the normal interaction forces, the influence of machining deviations on the thread load transfer law is clearly revealed. The generated axial load distribution feature sequence not only preserves the stress mutation signal induced by micro-defects but also provides high signal-to-noise ratio input data for subsequent accurate identification of fatigue hazard points through gradient convolution.

[0073] Further, the step of performing gradient convolution operation on the load distribution feature sequence to calculate the peak value of local field variables specifically includes: constructing a Gaussian first-order derivative convolution kernel, wherein the Gaussian first-order derivative convolution kernel is a discrete sequence of derivatives of a Gaussian function, used to extract the gradient features of the signal and suppress high-frequency noise; performing a sliding convolution operation on the load distribution feature sequence and the Gaussian first-order derivative convolution kernel in the time domain to output a gradient sequence reflecting the rate of change of the load distribution; searching for the magnitude maxima in the gradient sequence, locking the node position corresponding to the magnitude maxima in the time domain, and extracting the amplitude of the load distribution feature sequence at the node position as the peak value of the local field variables.

[0074] Specifically, a Gaussian first-order derivative convolution kernel is constructed. In this embodiment, for the non-stationary load signal induced by processing defects in bridge anchor bolts, the Gaussian first-order derivative convolution kernel refers to a discrete operator sequence generated by discretely sampling the first derivative of the Gaussian distribution function. The specific construction process is as follows: First, the window length L and smoothing standard deviation σ of the convolution kernel are preset. The value of the window length of the convolution kernel is determined based on the ratio of the nominal pitch of the bolt to the step size of the axial statistical interval, that is, the number of sampling points that can cover a complete pitch period is selected as the window length to ensure that the complete tooth profile features are included. The standard deviation is set to one-sixth of the window length. This value is set based on balancing the signal smoothing effect and the fidelity of the gradient features. Subsequently, the value of the Gaussian derivative function is calculated on the discrete integer time axis. The calculation logic is as follows: take the offset of the current sampling point relative to the center position, multiply it by the negative sign, divide it by the square of the standard deviation, and finally multiply it by the standard Gaussian distribution function value corresponding to the sampling point. The resulting discrete sequence is the Gaussian first derivative convolution kernel, which has a symmetrical peak and trough shape. It can effectively suppress discrete numerical noise caused by the non-uniformity of the finite element mesh in the load distribution sequence and highlight the abrupt gradient characteristics in the signal.

[0075] The load distribution feature sequence is subjected to a sliding convolution operation in the time domain with the Gaussian first derivative convolution kernel, outputting a gradient sequence reflecting the rate of change of the load distribution. The sliding convolution operation includes intermediate alignment: for a one-dimensional load distribution feature sequence of length 1000 rows and 1 column, a boundary copying and padding operation is first performed, that is, the value of the first data point is copied at the beginning of the sequence, and the value of the last data point is copied at the end, filling in a preset number of elements to avoid artificial gradient abrupt changes at signal edges, ensuring that the output dimension after convolution remains 1000 rows and 1 column. The specific execution logic of the sliding convolution operation is as follows: the Gaussian first derivative convolution kernel with a window length of 7 slides from the beginning of the load distribution feature sequence; at each sampling position, the 7 weight values ​​in the convolution kernel are element-wise multiplied with the 7 load values ​​in the corresponding window of the load sequence, and the multiplication results are summed to obtain the gradient prediction value corresponding to the current sampling point. The convolution kernel slides towards the end of the axis with a step size of one sampling point until the entire sequence is traversed. The resulting gradient sequence reflects the intensity of the load distribution at each axial position, with its positive and negative peaks representing the regions where the load rises or falls sharply along the bolt axis, respectively.

[0076] In a preferred embodiment, the step of performing gradient convolution operation on the load distribution feature sequence to calculate the peak value of local field variables includes: constructing a dual-scale Gaussian derivative convolution operator array containing a large-scale smoothing kernel and a small-scale feature kernel; using the large-scale smoothing kernel to perform macroscopic trend extraction on the load distribution feature sequence to generate a reference envelope sequence characterizing the axial force distribution of the threaded pair; using the small-scale feature kernel to extract high-frequency perturbation components in the load distribution feature sequence to generate a defect feature sequence; after trunculating the defect feature sequence with a significance threshold, performing element-level nonlinear weighted mapping with the reference envelope sequence to obtain an enhanced local gradient sequence, and extracting the peak value of local field variables from it.

[0077] Specifically, a dual-scale Gaussian derivative convolution operator array is constructed, comprising a large-scale smoothing kernel and a small-scale feature kernel. The process of constructing the large-scale smoothing kernel is as follows: based on the nominal pitch value of bridge bolts (e.g., 6 mm), the sampling span of the smoothing kernel is preset to 1.5 times the pitch (i.e., 9 mm). A long sequence operator with smooth transition is generated by calculating the first-order Gaussian derivative, used to suppress local mesh noise and lock the smooth attenuation trend of the load along the axial principal stress direction. The process of constructing the small-scale feature kernel is as follows: based on the resolution of the point cloud scan (e.g., 0.1 mm), the sampling span of the feature kernel is preset to 2 grid cell lengths (i.e., 0.4 mm). A spike-capturing operator with high sensitivity is generated by calculating the narrow-window Gaussian derivative.

[0078] A large-scale smoothing kernel is used to extract macroscopic trends from the load distribution feature sequence, generating a baseline envelope sequence characterizing the axial force distribution of the threaded joint. The computational logic involves performing a convolution and slip operation between the large-scale smoothing kernel and the normalized load distribution feature sequence. The baseline envelope sequence refers to the theoretical overall distribution curve reflecting the decreasing stress on each thread of a bridge bolt from the load end to the end. This sequence, through the large kernel smoothing effect, effectively filters out non-structural signal burrs caused by the micro-roughness of individual thread surfaces.

[0079] High-frequency perturbation components in the load distribution feature sequence are extracted using a small-scale feature kernel to generate a defect feature sequence. The computational logic is as follows: the small-scale feature kernel is convolved with the same load distribution feature sequence. Due to the extremely narrow operator window and its sensitivity to rate of change, this step can keenly capture minute load fluctuations caused by local geometric distortions mapped from discrete point clouds (such as non-uniform pitch transitions or local tooth flank collapse). The resulting defect feature sequence is a high-frequency oscillating vector containing a large number of local extrema.

[0080] A nonlinear weighted mapping is performed to obtain an enhanced local gradient sequence. The processing logic is as follows: First, the defect feature sequence is truncated based on a significance threshold, retaining only abrupt changes whose absolute value is greater than three times the sequence standard deviation. Then, the coordinates of the truncated feature points in the baseline envelope sequence are retrieved, and the envelope strength value at that location is extracted as a gain coefficient. The gain coefficient is then multiplied element-wise with the feature point value, resulting in a nonlinear amplitude amplification of the geometric defect signal located in a high-stress load-bearing region (such as the first thread). Through this composite mapping process, the "load background" and "morphological defects" are successfully mechanically coupled, generating a local gradient sequence that accurately indicates fatigue risk points.

[0081] By introducing a dual-scale convolution operator array and nonlinear weighted mapping, a high signal-to-noise ratio separation of micron-level machining defect signals from macroscopic load trends is achieved. This processing method effectively solves the problem that real defect signals in aspect ratio bolts are easily masked by the structural stress field, significantly enhancing the physical authenticity and identification sensitivity of fatigue weak point location.

[0082] The algorithm searches for the modulus maxima in the gradient sequence and locates the corresponding node position in the time domain. The modulus maxima refers to a feature point in the gradient sequence that has a local maximum absolute value and is representative in terms of energy. The specific search logic is as follows: Perform an absolute value operation on each value in the gradient sequence to generate a gradient modulus sequence. Iterate through the gradient modulus sequence and determine whether the value of the current sampling point is simultaneously greater than the values ​​of its preceding and following adjacent positions to identify local extrema. During this process, a significance threshold filtering based on machining tolerance constraints is performed. The preset logic for the significance threshold is as follows: Read the preset machining tolerance grade of the bridge bolts (e.g., 6g tolerance) and calculate the theoretical load gradient background value caused by the maximum allowable geometric fluctuation of this tolerance. Set the significance threshold to 1.5 times the mean of all values ​​in the gradient modulus sequence. This multiple is based on the fact that, in the case of large bridge anchorage, sensitivity tests on a sample set containing preset defects have shown that a threshold of 1.5 times the mean can effectively filter out minor burrs caused by finite element numerical rounding errors, while retaining the load abrupt changes caused by actual pitch steps or tooth flank distortion. By locking the maximum point that satisfies the above conditions, the node index position corresponding to the modulus maximum in the axial spatial domain can be accurately locked.

[0083] The amplitude of the load distribution feature sequence at the node location is extracted as the peak value of the local field variable. The peak value of the local field variable refers to the local stress control parameter that directly induces material fatigue damage. The specific extraction logic is as follows: using the index number of the modulus maxima sampling point locked in the preceding steps, the original one-dimensional load distribution feature sequence is accessed backtrackingly; the load intensity amplitude corresponding to the same index number in this sequence is extracted and defined as the peak value of the local field variable. This value comprehensively reflects the peak load intensity caused by microscopic geometric defects at this axial location and is a core variable for evaluating the fatigue life of bridge bolts.

[0084] By introducing a Gaussian first-derivative convolution kernel to extract gradients from the load distribution sequence, noise suppression and significant enhancement of local stress concentration characteristics are achieved. This processing method utilizes modulus maxima search logic to automatically locate the geometric points most severely affected by processing defects from complex axial load fluctuations. This not only eliminates the error of manually identifying fatigue weak points but also provides accurate physical load basis for subsequent fatigue life prediction based on damage evolution equations by extracting peak values ​​of local field variables, significantly improving the robustness of the prediction results.

[0085] Furthermore, the process of calculating the cumulative damage index based on the peak value of local field variables and the damage evolution equation, and outputting the predicted life of the target batch of high-strength bolts, includes: calling the mean stress correction model, substituting the peak value of the local field variables as the mean stress parameter into the mean stress correction model, and calculating the number of fatigue life cycles under the current load level; calculating the reciprocal of a single load cycle relative to the number of fatigue life cycles to obtain a single damage value; performing linear accumulation on the single damage values ​​under all load conditions throughout the entire life cycle to generate a cumulative damage index, and determining the total number of cycles when the cumulative damage index reaches a critical threshold, and using the total number of cycles as the predicted life.

[0086] Specifically, the mean stress correction model is invoked, and the peak value of the local field variable is substituted into the mean stress parameter to calculate the number of fatigue life cycles under the current load level. The mean stress correction model refers to a mathematical model used to evaluate the influence of asymmetric stress cycles on fatigue strength, such as the Morow correction model or the Gerber correction model. In this embodiment, the specific logic is as follows: First, the preset material basic fatigue parameters, including the fatigue strength coefficient and fatigue strength index, are retrieved. These parameters are pre-determined through standard fatigue tests on the same batch of alloy steel materials. Then, the mapping and analysis of stress parameters are performed: a spatial backtracking operation is performed: using the coordinate range of the axial statistical interval of the identified peak value of the local field variable, all contact nodes falling within this axial range are re-searched in the contact stiffness feature matrix; from this, the key node with the largest normal interaction value is selected, and the projected area of ​​the surface unit directly connected to the key node is extracted; the peak value of the local field variable is divided by the projected area corresponding to the key node, and the load dimension is converted into the stress intensity dimension to obtain a physical quantity characterizing the degree of local stress concentration at the micro-geometric defect, which is used as the maximum local stress response benchmark under the current working condition. Based on the cyclic characteristics (such as load ratio or peak-to-valley ratio) defined in the preset fatigue load spectrum, and combined with the maximum local stress response benchmark, the minimum local stress under the current working condition is calculated. Then, through arithmetic operations, the average stress component (i.e., the average of the maximum and minimum stresses) and stress amplitude (i.e., half of the difference between the maximum and minimum stresses) of the current cycle are obtained.

[0087] The specific method for calculating the fatigue life cycle count is as follows: Based on the selected mean stress correction model (taking the Moro model as an example), the fatigue strength coefficient of the material is corrected using the mean stress components (for example, mathematically, this is reflected as the fatigue strength coefficient minus the mean stress components, thereby reducing the effective fatigue resistance of the material). The corrected strength parameters, stress amplitude, and fatigue strength index are substituted into the standard stress-life equation, and the theoretical maximum number of cycles that the bolt can withstand before fatigue failure under the influence of the current geometric defects is obtained through exponentiation.

[0088] The single-cycle damage value is obtained by calculating the reciprocal of the number of fatigue life cycles relative to the single load cycle. The single-cycle damage value refers to the proportion of irreversible damage to the structural integrity of the bridge bolts caused by each stress cycle under the coupled action of a specific geometric defect and load. The specific calculation process is as follows: using a constant 1 as the dividend and the number of fatigue life cycles calculated in the previous steps as the divisor, perform arithmetic division. The quotient is the single-cycle damage value. The closer this value is to 1, the closer the transformed local contact stress caused by the current microscopic geometric defect is to the material's fatigue limit, and the greater the loss of bolt life due to a single load.

[0089] The single damage values ​​under all load conditions throughout the entire life cycle are linearly accumulated to generate a cumulative damage index. The total number of cycles required for the cumulative damage index to reach a critical threshold is determined, and this total number of cycles is used as the predicted life. The linear accumulation follows the Pamgren-Miner linear damage accumulation theory, which is a mechanical evaluation criterion that assumes that the damage to the structure under different stress levels is independent and linearly superimposed. The specific execution logic is as follows: First, a preset bridge service load spectrum is obtained, which defines the frequency of traffic load cycles of various amplitudes experienced by the bridge bolts during their service life. For each type of load condition in the load spectrum, the stress parameter conversion and single damage value calculation process are repeatedly performed using the aforementioned benchmark local contact stress combined with the load ratio of the current load condition. The single damage values ​​under all load conditions are weighted and accumulated according to their frequency of occurrence to generate a cumulative damage index. The critical threshold is a damage limit preset based on the engineering failure criterion, which is preset to 1.0 in this embodiment. The total number of load cycles required for the cumulative damage index to grow from zero to 1.0 is calculated, and this total number is defined as the predicted life. This process transforms the microscopic geometric defect characteristics obtained through optical detection into quantitative service life prediction indicators, providing data for the replacement cycle of bridge anchor bolts.

[0090] By incorporating the peak values ​​of local field variables reflecting processing defects into the mean stress correction model, this embodiment achieves cross-scale mechanical correlation from microscopic geometry to macroscopic fatigue damage. This processing method utilizes linear cumulative damage theory to couple complex random load spectra with real geometric defects, eliminating the life prediction bias caused by ideal geometry in traditional methods. The final predicted life accurately reflects the service safety differences of bolts in the target batch due to manufacturing process fluctuations, greatly improving the scientific rigor and accuracy of fatigue reliability assessment for critical bridge fasteners.

[0091] The fatigue life prediction method for high-strength bolted connections provided by this invention achieves full feature preservation of microscopic geometric defects in bridge anchor bolts in a numerical simulation environment by establishing a normal deviation mapping mechanism between discrete point cloud coordinate sequences and standard finite element mesh models. Using a finite element model reconstructed from a geometric defect mapping tensor, combined with a contact stiffness feature matrix calculated based on the penalty function method, it can accurately capture local contact stiffness singularities caused by cumulative pitch steps or tooth flank angle disturbances, solving the technical problem of failing to reflect the non-uniform distribution of contact stress under ideal geometric assumptions. By performing Gaussian first-derivative gradient convolution on the axial load distribution feature sequence and using modulus maxima search logic to eliminate numerical calculation noise, it achieves automated and accurate positioning of local stress concentration points induced by processing defects. Combining the average stress correction model and linear cumulative damage theory, this method successfully establishes a cross-scale quantitative correlation from physical measurement deviations to macroscopic service life, significantly improving the reliability and physical authenticity of life prediction for target batches of high-strength bolts under complex alternating load spectra, and providing a scientific numerical analysis basis for the safety assessment of large bridge fastening systems.

[0092] Example 2: This embodiment provides a fatigue life prediction system for high-strength bolted connections. This system is applied to the connection protection scenario of the anchorage zone at the cable end of a large suspension bridge, and performs fatigue service condition assessment on special high-strength alloy steel bolts. In this application scenario, because the anchor bolts are subjected to alternating loads and high stress levels for extended periods, micron-level processing deviations generated during manufacturing (such as waviness on the thread flank surface or non-uniform accumulation of the pitch) directly alter the load transfer path at the contact interface. This system, through the collaborative use of digital modules, transforms the unstructured morphological data of the bolt entity into local field variable characteristics that can drive fatigue prediction.

[0093] A fatigue life prediction system for high-strength bolted connections, referring to Figure 3The system includes: a data acquisition module for acquiring discrete point cloud coordinate sequences characterizing the physical morphology of a target batch of high-strength bolts; a geometric reconstruction module for reconstructing the discrete point cloud coordinate sequences using non-uniform rational B-spline interpolation to generate a continuous geometric topology function; a grid generation method for discretizing preset standard thread geometric design parameters to construct a standard finite element mesh model; a calculation of the difference in normal spatial position based on the continuous geometric topology function and the standard finite element mesh model, followed by matrix encapsulation to generate a geometric defect mapping tensor; and a node coordinate calculation based on the geometric defect mapping tensor of the standard finite element mesh model. The system is divided into four modules: a transformation module to construct a realistic finite element model; a simulation analysis module to construct a set of nonlinear numerical equilibrium equations based on the realistic finite element model using the penalty function method, and iteratively solve the normal interaction values ​​between each discrete contact node pair to generate a contact stiffness feature matrix; a load distribution feature sequence distributed along the axial direction based on the contact stiffness feature matrix, and gradient convolution operation to calculate the peak value of local field variables; and a life prediction module to calculate the cumulative damage index based on the peak value of local field variables and the damage evolution equation, and output the predicted life of the target batch of high-strength bolts.

[0094] Specifically, during the operation of bridge anchor bolts, the data acquisition module first reads the raw three-dimensional scattered point cloud data containing millions of spatial coordinates output by the optical measurement equipment deployed on the production line or testing station. Internally, this module constructs a three-dimensional KD tree spatial index structure, calculating the variance of each coordinate axis and performing spatial partitioning using the median to achieve hierarchical storage of massive amounts of data. The module traverses each data point, using the three-dimensional KD tree to search for a preset number (e.g., 20) of nearest neighbors and calculates the average Euclidean distance. By calculating the global mean and standard deviation, isolated measurement noise with a distance greater than the mean plus twice the standard deviation is removed. Finally, the module maps the cleaned point cloud from the Cartesian coordinate system to the cylindrical coordinate system, performs numerical ascending rearrangement according to the calculated global cumulative phase angle, and outputs a set of ordered discrete point cloud coordinate sequences that can describe the true helical morphology trajectory of the bridge bolts.

[0095] During execution, the geometric reconstruction module first uses the cumulative chord length parameterization method to assign a unique parameter node value to each point based on the physical chord length ratio between adjacent measurement points in the ordered sequence. This module then uses least-squares linear approximation logic to calculate the coordinates of the control vertices under a preset constraint of minimizing the sum of squared residuals. Combined with preset third-order basis functions, it constructs a continuous geometric topology function, which serves as the mathematical analytical benchmark describing the true tooth profile surface (including machining defects) of bridge bolts. Simultaneously, the module reads preset standard thread geometric design parameters (e.g., M100, 6mm pitch) and uses a hexahedral sweep algorithm to construct a defect-free standard finite element mesh model. By constructing normal rays pointing to the true morphological surface at ideal nodes and using Newton's iterative root-finding algorithm to calculate the spatial intersection points, the module calculates the coordinate deviation between the ideal nodes and the true intersection points and encapsulates it as a geometric defect mapping tensor corresponding to the dimension and the total number of nodes. Finally, the module corrects the node coordinates of the standard mesh through matrix superposition operations, constructing a true morphological finite element model that bears microscopic geometric deviations while maintaining the element connection topology unchanged.

[0096] During execution, the simulation analysis module first defines the nut's internal thread as the principal surface and the bolt's external thread as the slave surface on the finite element model of the actual morphology. It then presets a penalty stiffness coefficient (e.g., on the order of 10^13 Newtons per cubic meter) based on the bridge's material properties. The module detects the normal penetration depth of the slave surface nodes to the principal surface, accumulates the penalty stiffness terms corresponding to the penetrating nodes into the global stiffness matrix, and constructs a set of nonlinear equations describing the bridge's preload equilibrium state. Through Newton-Raphson iterative calculation, it outputs the normal interaction values ​​(i.e., contact force values) in the convergent state and fills them into the corresponding index positions of the global sparse diagonal matrix. The module generates a contact stiffness feature matrix. Further, it discretizes the entire axial length of the bolt into several spatial statistical intervals. By retrieving the node indices falling within these intervals and performing integral accumulation and normalization, it generates a one-dimensional load distribution feature sequence reflecting the axial load fluctuation characteristics. Subsequently, the module performs sliding convolution on this sequence using a Gaussian first-order derivative convolution kernel with a preset window length of 7, outputting a load change rate gradient sequence. Through searching for the maximum modulus and its significance filtering logic, it accurately identifies the axial position most severely affected by machining deviations (such as pitch steps) and extracts the load amplitude at that position as the peak value of the local field variable.

[0097] During execution, the life prediction module invokes a preset average stress correction model (such as the Moro correction model), using the extracted local field variable peaks as the average stress input. Combined with preset material fatigue strength parameters, it calculates the fatigue life cycle count under the current bridge service load level. This module utilizes linear damage accumulation theory (Pamgren-Miner theory) to calculate the damage contribution corresponding to a single load cycle and performs linear weighted accumulation based on a preset bridge service cycle load spectrum. When the cumulative damage index reaches a preset critical threshold of 1.0, the module determines and outputs the corresponding total number of cycles, which serves as the predicted life of the target batch of high-strength bolts in the current bridge application scenario.

[0098] Through the digital coordination of four functional modules—data acquisition, geometric reconstruction, simulation analysis, and life prediction—this system achieves non-destructive mapping and accurate prediction of fatigue service life from measured point cloud data of bridge bolts. By quantifying microscopic geometric defects generated during the manufacturing stage into the coordinate displacements of finite element nodes, the system overcomes the technical limitation of traditional simulation systems in capturing local singularities of contact stiffness. The introduction of gradient convolution and modulus maxima identification logic enables the system to automatically identify fatigue weak zones caused by manufacturing errors from massive amounts of data, greatly improving the robustness and scientific rigor of life prediction for key bridge fasteners under complex load spectra.

[0099] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention shall be determined by the scope of protection defined in the claims.

Claims

1. A method of predicting the fatigue life of a high-strength bolted joint, characterized by, The application relates to a method for predicting the service life of high-strength bolts. The method comprises the following steps: acquiring a discrete point cloud coordinate sequence representing the entity topography of a target batch of high-strength bolts; reconstructing the discrete point cloud coordinate sequence by a non-uniform rational B-spline interpolation method to generate a continuous geometric topology function; discretizing a preset standard thread geometric design parameter by a meshing method to construct a standard finite element mesh model; calculating a normal space position difference value and performing matrix encapsulation based on the continuous geometric topology function and the standard finite element mesh model to generate a geometric defect mapping tensor; performing node coordinate transformation on the standard finite element mesh model based on the geometric defect mapping tensor to construct a real topography finite element model; constructing a nonlinear numerical balance equation set by a penalty function method based on the real topography finite element model, and iteratively solving the normal interaction value between each discrete contact node pair to generate a contact stiffness characteristic matrix; calculating an axial distribution load distribution characteristic sequence based on the contact stiffness characteristic matrix, and performing gradient convolution operation on the load distribution characteristic sequence to calculate a local field variable peak value; 2. The method of predicting the fatigue life of a high-strength bolted joint according to claim 1, characterized in that, calculating a cumulative damage index based on the local field variable peak value and a damage evolution equation, and outputting the predicted service life of the target batch of high-strength bolts. The process of acquiring the discrete point cloud coordinate sequence representing the entity topography of the target batch of high-strength bolts comprises the following steps:

3. The method of predicting the fatigue life of a high-strength bolted joint according to claim 1, characterized in that, reading original three-dimensional scattered point cloud data output by an optical screening device, and constructing a three-dimensional KD tree space index structure; traversing each data point in the original three-dimensional scattered point cloud data, searching for K neighborhood points of the data point by using the three-dimensional KD tree space index structure, and calculating the average Euclidean distance of the data point to all K neighborhood points; calculating the mean and standard deviation of the average Euclidean distances of all data points, and removing outlier noise points with an average Euclidean distance greater than the mean plus twice the standard deviation; mapping the remaining data points after removing the noise points from a Cartesian coordinate system to a cylindrical coordinate system, and rearranging the data points in order according to the increasing order of the phase angle values of the spiral lines to generate an ordered discrete point cloud coordinate sequence. The generation process of the continuous geometric topology function comprises the following steps:

4. The method of predicting the fatigue life of a high-strength bolted joint according to claim 1, characterized by, calculating the parameter node vector of each data point in the discrete point cloud coordinate sequence by using an accumulated chord length parameterization method; constructing a linear approximation equation set based on a least square method, substituting the discrete point cloud coordinate sequence as the target observation value into the linear approximation equation set, and solving the control vertex coordinates that minimize the equation set residual sum of squares; constructing a parameterized continuous geometric topology function by using the control vertex coordinates, the parameter node vector and a preset basis function order, and the continuous geometric topology function is an analytical expression for describing the geometric trajectory of the thread surface. The generation process of the geometric defect mapping tensor comprises the following steps: The standard thread geometry parameter is read to generate an ideal spiral surface, and a hexahedral element is used to discretize the ideal spiral surface in space to construct a standard finite element grid model including a node coordinate matrix and an element connection relationship; discrete nodes on a contact surface in the standard finite element grid model are traversed, and a ray equation is constructed along a unit normal vector of each discrete node; the ray equation and the continuous geometric topology function are solved by a numerical root algorithm to obtain spatial intersection coordinates of an analytical surface defined by the ray equation and the continuous geometric topology function; a Euclidean distance vector between the spatial intersection coordinates and reference coordinates of the discrete nodes is calculated, and the Euclidean distance vector is decomposed into three Cartesian coordinate components, which are written into a sparse matrix corresponding to the node index to generate a geometric defect mapping tensor.

5. The method of predicting the fatigue life of a high-strength bolted joint according to claim 1, characterized in that, The construction process of the real topography finite element model includes: An initial node coordinate matrix of the standard finite element grid model is extracted, and a coordinate deviation matrix stored in the geometric defect mapping tensor is extracted; a matrix linear superposition operation is performed on the initial node coordinate matrix and the coordinate deviation matrix to generate an updated node coordinate matrix; the spatial node definition of the grid model is reset by using the updated node coordinate matrix, the micro-geometric deviation is introduced on the premise of keeping the grid topology connection relationship unchanged, and a real topography finite element model is constructed.

6. The method of predicting the fatigue life of a high-strength bolted joint according to claim 1, characterized in that, The process of constructing a nonlinear numerical balance equation set by a penalty function method and iteratively solving the normal interaction value between each pair of discrete contact nodes includes: A contact major surface and a contact minor surface are defined in the real topography finite element model, and a penalty stiffness coefficient is set; the normal penetration depth of a node on the contact minor surface relative to the contact major surface is detected, and for a node pair with a normal penetration depth greater than zero, the penalty stiffness coefficient is added to the corresponding diagonal element position of the overall stiffness matrix; the nonlinear numerical balance equation set including the penalty stiffness coefficient is constructed, and the Newton-Raphson algorithm is used to iteratively solve the nonlinear numerical balance equation set until the norm of the residual vector meets the preset convergence criterion; the product of the penalty stiffness coefficient and the normal penetration depth in the convergence state is output as the normal interaction value.

7. The method of predicting the fatigue life of a high-strength bolted joint according to claim 1, characterized by, The process of calculating an axial distribution load distribution feature sequence based on the contact stiffness feature matrix includes: A global sparse diagonal matrix with a dimension equal to the total number of nodes of the real topography finite element model is constructed; all contact node pairs are traversed, the normal interaction value of each contact node pair is taken as a feature element, and the feature element is filled into the diagonal line position corresponding to the node index in the global sparse diagonal matrix to generate a contact stiffness feature matrix; non-zero diagonal elements in the contact stiffness feature matrix are extracted, and the axial length of the thread is discretized into a plurality of continuous axial statistical intervals; all node indexes falling into each of the axial statistical intervals are retrieved, the normal interaction values corresponding to the node indexes in the contact stiffness feature matrix are extracted, and the extracted normal interaction values are integrated and summed to generate a one-dimensional load distribution feature sequence.

8. The method of predicting the fatigue life of a high-strength bolted joint according to claim 1, characterized by, The gradient convolution operation is performed on the load distribution characteristic sequence, and a local field variable peak value is calculated, specifically including: A Gaussian first derivative convolution kernel is constructed, which is a derivative discrete sequence of a Gaussian function, used for extracting gradient features of the signal and suppressing high-frequency noise; a sliding convolution operation is performed on the load distribution characteristic sequence and the Gaussian first derivative convolution kernel in the time domain, and a gradient sequence reflecting the rate of change of the load distribution is output; a modulus maximum value in the gradient sequence is searched, a node position corresponding to the modulus maximum value in the time domain is locked, and an amplitude value of the load distribution characteristic sequence at the node position is extracted as a local field variable peak value.

9. The method of predicting the fatigue life of a high-strength bolted joint according to claim 1, characterized by, The process of calculating a cumulative damage index based on the local field variable peak value and combining a damage evolution equation, and outputting a predicted life of the target batch of high-strength bolts includes: The average stress correction model is called, the local field variable peak value is taken as an average stress parameter, the fatigue life cycle number under the current load level is calculated, the reciprocal of a single load cycle with respect to the fatigue life cycle number is calculated to obtain a single damage value, and the single damage values under all load conditions in the whole life cycle are linearly accumulated to generate a cumulative damage index, and the total cycle number when the cumulative damage index reaches a critical threshold is determined as the predicted life.

10. A high-strength bolted joint fatigue life prediction system, characterized by, It includes: A data acquisition module is configured to acquire a discrete point cloud coordinate sequence representing an entity morphology of the target batch of high-strength bolts; A geometry reconstruction module is configured to reconstruct the discrete point cloud coordinate sequence by a non-uniform rational B-spline interpolation method to generate a continuous geometric topology function, and to discretize a preset standard thread geometry design parameter by a meshing method to construct a standard finite element mesh model; Based on the continuous geometric topology function and the standard finite element mesh model, a normal spatial position difference value is calculated and matrix-encapsulated to generate a geometric defect mapping tensor; Based on the geometric defect mapping tensor, a node coordinate transformation is performed on the standard finite element mesh model to construct a real morphology finite element model; A simulation analysis module is configured to construct a nonlinear numerical balance equation set by a penalty function method based on the real morphology finite element model, and to iteratively solve a normal interaction value between each discrete contact node pair to generate a contact stiffness feature matrix, calculate a load distribution characteristic sequence distributed along an axial direction based on the contact stiffness feature matrix, and perform a gradient convolution operation on the load distribution characteristic sequence to calculate a local field variable peak value. A life prediction module is configured to calculate a cumulative damage index based on the local field variable peak value and combining a damage evolution equation, and output a predicted life of the target batch of high-strength bolts.

Images