Asphalt pavement compression analysis system based on dynamic simulation

Abstract

The present application relates to the technical field of computer aided engineering, specifically to an asphalt pavement compression resistance analysis system based on dynamic simulation, which comprises a tire load discretization module, which is used to generate a contact mechanics pre-calculation and map a ground pressure distribution cloud map according to three-dimensional scanning geometric data of a tire, extract node pressure values in the ground pressure distribution cloud map, and arrange the node pressure values according to tire pattern geometric characteristics.In the present application, the cumulative degree of internal damage of materials can be characterized from the dual dimensions of energy dissipation and shape distortion, the growth rate slope of the hysteresis loop geometric distortion index relative to the number of loading cycles can be monitored, and the slope mutation coordinates can be located, thereby achieving sensitive capture of the compression fatigue life ending inflection point.This discrimination logic based on evolution rate mutation is more accurate than simple stress threshold determination, thereby improving the reliability of asphalt pavement structure compression fatigue life prediction.

Description

Technical Field

[0001] This invention relates to the field of computer-aided engineering technology, and in particular to an asphalt pavement compressive strength analysis system based on dynamic simulation. Background Technology

[0002] Computer-aided engineering technology utilizes computer software and powerful computing capabilities to simulate, analyze, and optimize engineering problems, aiming to solve various physical and mechanical challenges from product design to manufacturing processes.

[0003] Current computer-aided engineering techniques, when simulating road structures, typically simplify the tire-road contact as a uniformly distributed circular or rectangular static load, neglecting the local stress concentration effects at the edges of tire tread blocks and the changes in contact pressure gradients. This results in a flattened stress distribution field in the simulated road surface. Furthermore, during dynamic response analysis, moving loads crossing mesh boundaries can easily induce spurious high-frequency noise, masking the true mechanical response signal and affecting the reliability of the calculation results. Simultaneously, traditional fatigue life assessment methods often rely on empirical formulas based on maximum tensile strain or peak stress, focusing only on extreme responses at specific moments while ignoring the continuous evolution of the viscoelastic hysteresis characteristics of materials under cyclic loading. Simply relying on peak values ​​makes it difficult to effectively distinguish between recoverable viscoelastic deformation and irreversible damage accumulation. Therefore, improvements are needed. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and propose an asphalt pavement compressive strength analysis system based on dynamic simulation.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: An asphalt pavement compressive strength analysis system based on dynamic simulation includes: The tire load discretization module is used to perform contact mechanics pre-calculation based on the three-dimensional scanning geometric data of the tire and map it to generate a ground pressure distribution cloud map. It extracts the node pressure values ​​in the ground pressure distribution cloud map, arranges the node pressure values ​​according to the geometric features of the tire tread pattern, and generates a high-resolution pressure matrix. The dynamic pressure field loading module is used to distribute pressure values ​​to the road surface finite element mesh nodes according to the high-resolution pressure matrix to generate nodal force load vectors, calculate the load transfer ratio between adjacent mesh nodes according to the nodal force load vectors, and generate a smooth dynamic load sequence. The viscoelastic hysteresis loop extraction module is used to perform time-step integral calculations on the viscoelastic constitutive parameters of asphalt pavement according to the smooth dynamic load sequence, draw and generate a closed hysteresis loop curve, and perform least squares fitting on the closed hysteresis loop curve in combination with the current dynamic modulus and phase angle parameters to generate an ideal viscoelastic elliptical geometry. The compressive fatigue life prediction module is used to calculate the Frescher distance difference value on the geometric topology based on the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry, obtain the hysteresis loop geometric distortion index, calculate the slope of the growth rate of the hysteresis loop geometric distortion index relative to the number of loading cycles, locate the coordinates of the slope change, and generate the inflection point of the end of compressive fatigue life.

[0006] Preferably, the step of obtaining the high-resolution pressure matrix is ​​as follows: Based on the three-dimensional scanning geometric data of the tire, the boundaries of the tread blocks and polygonal patches are analyzed, the scanning coordinate axis direction is corrected, the patch with the normal vector facing the ground is located, the contact units are aggregated according to the ground projection, the pressure values ​​are distributed according to the contact ratio of the patch, and a ground pressure distribution cloud map is generated. Based on the ground pressure distribution cloud map, the minimum and maximum values ​​of the color level mapping table are read, and the color level value of each pixel is converted into pressure units according to the corresponding index from the pixel center coordinate to the ground coordinate. Non-contact pixels are eliminated and written in the order of grid number to generate node pressure values. Based on the node pressure values, an index table is constructed according to the block numbering order of the tire tread geometry and the groove centerline sequence. First, it is sorted radially from the center of the tire crown to the tire shoulder, and then sorted circumferentially clockwise. The node pressure values ​​of each row and column are filled in sequentially to generate a high-resolution pressure matrix.

[0007] Preferably, the step of obtaining the nodal force load vector is as follows: Based on the high-resolution pressure matrix, the vehicle's center of gravity position and roll angle parameters are read, and a projection index from the tire contact point to the road surface finite element mesh node is established. The pressure values ​​are then distributed to the road surface finite element mesh node using bilinear interpolation mapping to generate the nodal force load vector.

[0008] Preferably, the step of obtaining the smooth dynamic load sequence is as follows: Based on the nodal force load vector, read the adjacent grid node pairs according to the grid topology, calculate the pressure difference of each pair of nodes and the ratio of the total pressure as a local proportion, write it according to the grid edge number and unify the range to generate the load transfer proportion; Based on the load transfer ratio, the nodal force load vector is weighted and distributed among adjacent grid nodes according to the time step index. The load transfer ratio is used to smooth the abrupt changes between adjacent time steps and rearrange the node order. The nodal force load vectors of each time step are then spliced ​​together to generate a smooth dynamic load sequence.

[0009] Preferably, the step of obtaining the closed hysteresis loop curve is as follows: Based on the smooth dynamic load sequence, the time step and start and end time of the viscoelastic constitutive parameters of the asphalt pavement are read. The load increment of each time step is calculated according to the time index and the stress increment and strain increment are updated synchronously. The stress and strain values ​​are recorded in the order of the integration point number and the cycle number is marked to generate a sequence of stress and strain values ​​at key integration points. Based on the stress and strain numerical sequence of the key integration points, the start and end points of each cycle are located according to the cycle number and the corresponding intervals are extracted. The stress and strain values ​​are paired one by one according to the same time index and the curves are connected in chronological order. Discrete points outside the cycle are eliminated and the beginning and end are connected to meet the closure condition, thus generating a closed hysteresis loop curve.

[0010] Preferably, the steps for obtaining the ideal viscoelastic elliptic geometry are as follows: Based on the closed hysteresis loop curve, the dynamic modulus parameter and phase angle parameter are called to set the elliptical constraint coordinate system and initialize the major axis direction and minor axis direction. The deviation to the constraint coordinate system is calculated segment by segment according to the curve and the segment weight is set according to the dynamic modulus parameter. The major axis length, minor axis length and rotation angle are updated round by round until the total deviation no longer decreases, and an ideal viscoelastic elliptical geometry is generated.

[0011] Preferably, the step of obtaining the hysteresis loop geometric distortion index is as follows: Based on the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry, the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry are synchronously resampled at arc length intervals to establish a matching path table. The coordinates of the corresponding sampling points are matched step by step and the Euclidean distance between each pair of matching points is calculated. The minimum value of the maximum point distance is extracted along the path to generate the Fraser distance difference value. The hysteresis loop geometric distortion index is calculated based on the Fraser distance difference value.

[0012] Preferably, the step of obtaining the inflection point of the end of the compressive fatigue life is as follows: Based on the hysteresis loop geometric distortion index, the hysteresis loop geometric distortion index of consecutive cycles is differentially divided according to the number of loading cycles to obtain the growth rate slope sequence of the hysteresis loop geometric distortion index. When the growth rate of the hysteresis loop geometric distortion index first exceeds the preset damage threshold, and the rate remains above the preset damage threshold in the next two consecutive cycles, the coordinate point corresponding to the cycle is determined as the inflection point of the end of the compressive fatigue life.

[0013] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, contact mechanics pre-calculation is performed using three-dimensional scanning geometric data of the tire, and a ground pressure distribution cloud map is generated. This captures the non-uniform contact characteristics caused by the tire tread structure. The nodal pressure values ​​are arranged according to the tread geometry to generate a high-resolution pressure matrix, ensuring that the simulated boundary conditions closely match the actual road surface stress state. Bilinear interpolation is used to distribute the pressure values ​​to the road surface finite element mesh nodes and calculate the load transfer ratio between adjacent nodes, generating a smooth dynamic load sequence. This processing method eliminates numerical oscillations and stress singularities that may be caused by discretized meshes during dynamic loading, ensuring the convergence and numerical stability of the dynamic calculations. A closed loop is plotted by performing time-step integral operations on the viscoelastic constitutive parameters of the asphalt pavement. Hysteresis loop curves are fitted with dynamic modulus and phase angle to generate an ideal viscoelastic elliptical geometry, establishing an intuitive benchmark for comparing the ideal and damaged states of the material. The Fréchet distance difference between the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry on the geometric topology is calculated, and the hysteresis loop geometric distortion index is quantified. This index can characterize the accumulation of internal damage in the material from both energy dissipation and shape distortion dimensions. By monitoring the growth rate slope of the hysteresis loop geometric distortion index relative to the number of loading cycles and locating the coordinates of slope abrupt changes, the inflection point of the end of compressive fatigue life is accurately captured. This discrimination logic based on abrupt changes in evolution rate is more accurate than simple stress threshold determination, thereby improving the reliability of predicting the compressive fatigue life of asphalt pavement structures. Attached Figure Description

[0014] Figure 1 This is a system flowchart of the present invention. Detailed Implementation

[0015] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0016] Please see Figure 1 The present invention provides a technical solution: an asphalt pavement compressive strength analysis system based on dynamic simulation, comprising: The tire load discretization module is used to perform contact mechanics pre-calculation based on the three-dimensional scanning geometric data of the tire and map it to generate a ground pressure distribution cloud map. It extracts the nodal pressure values ​​in the ground pressure distribution cloud map, arranges the nodal pressure values ​​according to the geometric features of the tire tread pattern, and generates a high-resolution pressure matrix. The dynamic pressure field loading module is used to distribute pressure values ​​to the road surface finite element mesh nodes according to the high-resolution pressure matrix to generate nodal force load vectors, calculate the load transfer ratio between adjacent mesh nodes based on the nodal force load vectors, and generate a smooth dynamic load sequence. The viscoelastic hysteresis loop extraction module is used to perform time-step integral calculations on the viscoelastic constitutive parameters of asphalt pavement according to a smooth dynamic load sequence, draw and generate a closed hysteresis loop curve, and perform least squares fitting on the closed hysteresis loop curve in combination with the current dynamic modulus and phase angle parameters to generate an ideal viscoelastic elliptical geometry. The compressive fatigue life prediction module is used to calculate the Frescher distance difference value on the geometric topology based on the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry, obtain the hysteresis loop geometric distortion index, calculate the slope of the growth rate of the hysteresis loop geometric distortion index relative to the number of loading cycles, locate the coordinates of the slope change, and generate the inflection point of the end of the compressive fatigue life.

[0017] The steps for obtaining a high-resolution pressure matrix are as follows: Based on the three-dimensional scanning geometric data of the tire, the boundaries of the tread blocks and polygonal patches are analyzed, the scanning coordinate axis direction is corrected, the patch with the normal vector facing the ground is located, the contact units are aggregated according to the ground projection, the pressure values ​​are distributed according to the contact ratio of the patch, and a ground pressure distribution cloud map is generated. Based on the ground pressure distribution cloud map, the minimum and maximum values ​​of the color level mapping table are read. According to the corresponding index from the pixel center coordinate to the ground coordinate, the color level value of each pixel is converted into pressure unit. Non-contact pixels are eliminated and written in the order of grid number to generate node pressure values. Based on the node pressure values, an index table is constructed according to the block numbering order of the tire tread geometry and the groove centerline sequence. First, it is sorted radially from the center of the tire crown to the tire shoulder, and then sorted circumferentially clockwise. The node pressure values ​​of each row and column are filled in sequentially to generate a high-resolution pressure matrix.

[0018] Specifically, based on the 3D scanning geometric data of the tire, a point cloud file or mesh file containing all discrete points on the tire surface is read. The vertex list and face index list in the file are traversed, and the coordinates of the three vertices of each polygon face are extracted. The coordinates of the geometric center of each face are calculated as the face centroid. Principal component analysis is used to calculate the covariance matrix of the overall point cloud. Eigenvalue decomposition is performed on the covariance matrix to obtain the eigenvector corresponding to the largest eigenvalue. The direction of this eigenvector is defined as the rotation axis direction of the tire. The plane perpendicular to the rotation axis and passing through the geometric center of the tire is set as the equatorial plane. The origin of the coordinate system is translated to the intersection of the equatorial plane and the rotation axis. The Z-axis is adjusted to be perpendicular to the road surface direction. The normal vector of all faces is calculated, and a normal vector filtering threshold is set. For example, the value is This value is set empirically, representing the condition that the cosine of the angle between the normal vector and the direction of gravity must be less than this threshold to be considered as facing the ground. A set of patches meeting this condition is selected, and a two-dimensional projection mesh of the road surface contact area is created, with the mesh resolution set to [value missing]. The selected facets are projected along the Z-axis onto a two-dimensional projection grid plane. The projected area of ​​each facet within the grid is calculated. A polygon clipping algorithm is used to calculate the area of ​​the overlapping region between the facet projection and the grid cells. Finally, the facet contact ratio coefficient is calculated. The calculation formula is: in, The surface contact ratio coefficient represents the weight of the surface's contribution to the mesh pressure. The area of ​​the overlapping region between the patch projection and the mesh element is calculated using the shoelace formula from the vertex coordinates of the clipped polygon. The total area of ​​the grid cells, i.e. Obtain the total load per wheel set by the vehicle. For example, set to The total load is distributed to each contact unit according to the contact ratio coefficient, and the unit pressure value is calculated. The calculation formula is: in, The pressure value assigned to the current grid cell. For the total load of a single wheel, The calculated pressure value is mapped to a preset color space by summing the contact ratio coefficients of all effective contact surfaces, for example, by using linear interpolation from blue to red, to generate a color image containing pressure information and generate a ground pressure distribution cloud map.

[0019] Based on the ground pressure distribution cloud map, the header information of the image file is analyzed to obtain the width and height pixel count of the image, the color legend area in the image is located, and the lowest pressure value corresponding to the color scale bar in the legend is identified. With the highest pressure value For example, setting for , for A linear mapping function between RGB color values ​​and pressure values ​​is established. Then, each pixel in the ground pressure distribution cloud map is traversed to extract the red, green, and blue channel values ​​of the current pixel. Calculate the relative position ratio of the current pixel color in the color scale bar. The physical pressure value represented by the pixel is calculated based on the relative position ratio. The calculation formula is: in, This is the pressure unit value converted from the current pixel value. This represents the minimum boundary of the color level mapping table. This represents the normalized position of the pixel color value in the color gradient, with a value ranging from 0 to 1. Set a non-contact noise filtering threshold for the maximum value boundary of the color level mapping table. The threshold is set to the maximum pressure value. ,Right now , calculate and Perform a comparison, if Less than These are identified as background noise or non-contact areas and are removed, while valid contact pixels are retained, based on their row and column coordinates in the image. And the scale between image resolution and actual physical size (For example ), calculate the ground physical coordinates corresponding to the pixel center point Find the finite element mesh node number that is closest to the physical coordinates of the ground. The converted pressure values ​​are assigned to the corresponding grid nodes. All nodes containing effective pressure values ​​are arranged in ascending order of grid number and organized into a one-dimensional array to generate node pressure values.

[0020] Based on the node pressure values, the geometric parameters of the tread grooves in the tire design drawings are read, the centerline coordinates of the longitudinal main grooves are identified, the tire contact patch area is divided into several independent tread block areas, and each block is assigned a unique number. All nodes that record pressure values ​​are traversed, and the polar diameter of each node relative to the tire contact patch center point is calculated. and polar angle Based on the relative position of the node coordinates to the centerline of the groove, each node is classified into its corresponding tread block number. A node index table containing block number, polar diameter, polar angle, and pressure value is constructed. A two-level sorting operation is performed. First, a primary sorting is performed, using the tire radial direction as a reference, and then the nodes in the index table are sorted according to their polar diameter. From this point on the center of the fetal crown ( (towards the shoulder) (Increase) Sort in ascending order. If nodes are located in the same radial position or the difference in polar diameter is within the set tolerance, Within the range, a secondary sorting is performed, based on the tire circumference and according to the polar angle. The nodes are sorted in ascending order clockwise. A two-dimensional matrix container is created based on the sorted node order. The number of rows in the matrix corresponds to the number of radial discrete points, and the number of columns corresponds to the number of circumferential discrete points. The sorted node pressure values ​​are then filled into the matrix cells sequentially. For non-contact areas caused by tread grooves, zero values ​​are filled into the corresponding positions in the matrix. The matrix edges are smoothed to eliminate the jagged effect caused by discretization. Finally, a numerical matrix that can reflect the fine structure of the tire tread is constructed, generating a high-resolution pressure matrix.

[0021] The steps for obtaining the nodal force load vector are as follows: Based on the high-resolution pressure matrix, the vehicle's center of gravity position and roll angle parameters are read, and a projection index from the tire contact point to the road surface finite element mesh node is established. The pressure values ​​are then distributed to the road surface finite element mesh node using bilinear interpolation mapping to generate the nodal force load vector.

[0022] Specifically, based on the high-resolution pressure matrix, the vehicle's center of gravity position parameter is set as the center of gravity height. Roll angle parameters The single-wheel load correction factor caused by roll is calculated based on the load transfer principle in vehicle dynamics. The calculation formula is: in, This is a load correction factor used to adjust static loads. The roll stiffness of the suspension system is set to [value]. This value is selected based on common heavy-duty truck suspension parameters. The vehicle roll angle (in radians). Set the vehicle track width to 1. , The curb weight of the vehicle is set as , Let be the acceleration due to gravity, and take . This coefficient is used to perform a global multiplicative scaling update on all values ​​in the high-resolution pressure matrix, establishing a global coordinate system for the pavement finite element mesh. Each pavement finite element node is then traversed to obtain its planar coordinates. Inversely transform it to the local coordinate system of the high-resolution pressure matrix. Determine the indices of the four corner points of the pressure matrix element into which the coordinate falls, respectively. , , and Read the pressure values ​​at these four corner points. Calculate interpolation weights and Perform bilinear interpolation to calculate nodal pressure. The calculation formula is: in, This represents the pressure value mapped to the road surface node. The normalized distance weights of the sampling points in the horizontal direction relative to the left grid line. The normalized distance weights of the sampling points in the vertical direction relative to the grid lines below. Given the known pressure values ​​at the four corner points, the calculated pressure value is multiplied by the control area corresponding to that finite element node. The vertical force value of the node is obtained, and the force values ​​of all nodes are arranged in order of node number to generate the node force load vector.

[0023] The steps for obtaining a smoothed dynamic load sequence are as follows: Based on the nodal force load vector, read the adjacent grid node pairs according to the grid topology, calculate the pressure difference of each pair of nodes and the ratio of the total pressure as a local proportion, write it according to the grid edge number and unify the range to generate the load transfer proportion; Based on the load transfer ratio, the nodal force load vector is weighted and distributed among adjacent grid nodes according to the time step index. The load transfer ratio is used to smooth the abrupt changes between adjacent time steps and rearrange the node order. The nodal force load vectors of each time step are then spliced ​​together to generate a smooth dynamic load sequence.

[0024] Specifically, based on the nodal force load vector, the node connection table of the pavement finite element model is read, all mesh edges consisting of two nodes are identified, an adjacency matrix is ​​established to store the indices of adjacent node pairs, and each mesh edge is traversed to obtain the node numbers of the two nodes it connects to. and Extract the load values ​​of these two nodes from the nodal force load vector obtained in the previous steps. and Set the minimum effective load threshold For example, set to This value is set based on the sensor's noise floor level. If the load transfer ratio of the node pair is not specified, then set it to 0 to avoid division by zero error; otherwise, calculate the load transfer ratio. The calculation formula is: in, A dimensionless scaling factor to reflect the load gradient variation between adjacent nodes. and These are the vertical load force values ​​for the two nodes in a pair of adjacent nodes, with the following symbols: This indicates the absolute value operation; after calculation, it applies to all edges. Perform a normalization check; if any values ​​exceed the preset range... Outliers are truncated to boundary values. The calculation results are serialized and stored according to the unique identifier ID of the mesh edge. A lookup table containing edge ID, node A number, node B number and transfer ratio is constructed to generate the load transfer ratio.

[0025] Set the time step for dynamic loading based on the load transfer ratio. Total simulation time Initialize time step index For each time step, the simulated tire's position on the road surface is determined, and the nodal force load vector is updated based on the contact area after movement to obtain the initial load vector at the current moment. To eliminate the abrupt loading caused by mesh discretization, the aforementioned load transfer ratio is utilized. Smooth the adjacent nodes and traverse all adjacent node pairs. Calculate the smoothing correction amount based on its corresponding transmission ratio. The calculation formula is: in, This refers to the load correction amount that needs to be balanced between nodes. As a smoothing factor, set to This value was determined through preliminary trials to ensure the stability of the numerical evolution. This represents the load transfer ratio for this node pair. and For the current time step node and The original load values ​​are updated to the node loads. and After completing the smooth iteration of all node pairs, the node load values ​​of the current time step are arranged in global node number order to form a column vector. The column vectors of all time steps are then horizontally concatenated in chronological order to construct a dimension-1. The dynamic load matrix, where This represents the total number of road surface nodes. Given the total number of time steps, a smooth dynamic load sequence is generated.

[0026] The steps to obtain a closed hysteresis loop curve are as follows: Based on the smooth dynamic load sequence, the time step and start and end time of the viscoelastic constitutive parameters of the asphalt pavement are read. The load increment of each time step is calculated according to the time index and the stress increment and strain increment are updated synchronously. The stress and strain values ​​are recorded in the order of integration point number and the cycle number is marked to generate the stress and strain value sequence of key integration points. Based on the stress and strain numerical sequence of key integration points, the start and end points of each cycle are located according to the cycle number and the corresponding intervals are extracted. Stress and strain values ​​are paired one by one according to the same time index and the curves are connected in chronological order. Discrete points outside the cycle are removed and the beginning and end are connected to meet the closure condition, thus generating a closed hysteresis loop curve.

[0027] Specifically, the time step for reading the viscoelastic constitutive parameters of asphalt pavement is based on the smoothed dynamic load sequence. With start and end times, for example, setting the time step to... This value depends on the loading frequency. The sampling theorem is set to ensure that each cycle contains at least 50 sampling points, the start and end times are set to cover the complete vehicle passage process, the historical variables for stress and strain calculation are initialized, each time node in the smoothed dynamic load sequence is traversed, and the load value at the current moment is extracted. Combined with the effective contact area of ​​the road surface stress zone Calculate the vertical compressive stress at the current moment. The corresponding viscoelastic strain is calculated using a recursive integration algorithm based on the generalized Maxwell model or the generalized Kelvin model. The calculation formula is: in, The strain value is calculated at the current time step. The glassy compliance of asphalt mixtures represents the instantaneous elastic response. This is the number of terms in the Prony series, typically ranging from 5 to 9 to cover wide frequency range characteristics. For the first The hysteresis time constant of the term, For the first The compliance coefficient of the term, For the first The term refers to the historical integral variable at the previous time step, used to record the memory effect of the material. For the input stress at the current time step, synchronously update the historical variables of all Prony series terms. For use in the next step of calculation, based on the current accumulated time. With loading cycle The relationship is used to calculate the cycle number to which the current data point belongs. The calculation formula is: , where the symbol This indicates rounding down. The calculated stress values, strain values, timestamps, and cycle numbers are organized by column and stored in a pre-allocated dynamic array. For key locations in the pavement structure that are prone to fatigue cracking, such as the bottom of the surface layer or the center of the wheel track, the corresponding integration point data is extracted to generate a sequence of stress and strain values ​​at key integration points.

[0028] Based on the stress and strain numerical sequence at key integration points, the cycle number column in the sequence is traversed to identify the index positions where the sequence number jumps, thereby locating the starting index of each loading cycle. and end index Extract the stress data arrays from the complete cycle one by one. and strain data array Check the data integrity for each loop; if the number of data points is less than the preset minimum number of points threshold... (For example, if the value is set to 20 points), then it is considered an invalid loop and is discarded. The extracted valid loop data is then subjected to closure detection, and the Euclidean distance between the first and last data points is calculated. The calculation formula is: in, and These are the stress and strain values ​​at the end of the cycle. and Set a closure tolerance threshold for the stress and strain values ​​at the start of the cycle. The threshold value is taken as the current cyclic stress amplitude. ,like This indicates the presence of drift or residual deformation. A linear drift correction is performed to force closure. The correction formula is as follows: in, and For the revised first Stress and strain values ​​at each point and The original measured values ​​are subtracted to achieve overlap by subtracting the deviation that accumulates linearly over time. The corrected stress and strain data are then paired one by one according to the time series, and discrete points far from the cluster center caused by noise are removed. The processed data points are then connected in chronological order to form a closed polygonal path, generating a closed hysteresis loop curve.

[0029] The steps to obtain an ideal viscoelastic elliptic geometry are as follows: Based on the closed hysteresis loop curve, the dynamic modulus parameters and phase angle parameters are called to set the elliptical constraint coordinate system and initialize the major axis and minor axis directions. The deviation from the constraint coordinate system is calculated segment by segment according to the curve and the segment weight is set according to the dynamic modulus parameters. The major axis length, minor axis length and rotation angle are updated round by round until the total deviation no longer decreases, thus generating an ideal viscoelastic elliptical geometry.

[0030] Specifically, based on the closed hysteresis loop curve, the dynamic modulus of the asphalt pavement at the current temperature and loading frequency is read. and phase angle Using these two parameters, an initial ideal elliptical model is constructed, and the parametric equations of the ellipse are defined based on the stress amplitude. Initialize the major axis by taking half of the maximum stress value in the closed curve. and short axis Initial settings and A local coordinate system is established with the stress-strain geometric center as the origin. The objective function is defined as the sum of squared normal distances from each point on the closed hysteresis loop curve to the ideal elliptical trajectory. The Levenberg-Marquardt iterative algorithm is used to fine-tune the geometric parameters of the ellipse. During the iteration process, the geometric parameters of each sampling point on the closed hysteresis loop are calculated. Algebraic distance deviation to the current elliptic model The calculation formula is: in, This represents the total residual value. The total number of sampling points. Let be the rotation angle of the principal axis of the ellipse relative to the stress axis, initialized to half of the phase angle. and For the semi-major axis and semi-minor axis lengths to be optimized, and The coordinates of the sampled points after centralization are set, and the iteration termination condition is set to the residual change rate being less than 1 / 2 for two consecutive iterations. This threshold is set based on the convergence accuracy requirements of the numerical calculation and is continuously updated through iteration. , and The values ​​are calculated until the theoretical ellipse that best matches the actual hysteresis loop geometry is found. The ellipse parameters at the final convergence are recorded to generate the ideal viscoelastic elliptical geometry.

[0031] The steps to obtain the geometric distortion index of the hysteresis loop are as follows: Based on the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry, the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry are synchronously resampled at arc length intervals to establish a matching path table. The coordinates of the corresponding sampling points are matched step by step and the Euclidean distance between each pair of matching points is calculated. The minimum value of the maximum point distance is extracted along the path to generate the Fraser distance difference value. Based on the Frescher distance difference value, the geometric distortion index of the hysteresis loop is calculated using the following formula: ; in, The geometric distortion index of the hysteresis loop. Let F be the Fraser distance difference value, representing the minimum and maximum topological distances between the closed hysteresis loop curve and the ideal viscoelastic elliptic geometry. Let be the arc length of an ideal viscoelastic elliptic geometry, obtained by integrating over the circumferential coordinates of the ideal ellipse. The area enclosed by the closed hysteresis loop curve is obtained by line integration in the stress-strain coordinate plane. The area enclosed by an ideal viscoelastic elliptic geometry is calculated by multiplying the major and minor axes of the ideal ellipse.

[0032] Specifically, based on the closed hysteresis loop curve and the ideal viscoelastic elliptic geometry, the original stress-strain data points of the closed hysteresis loop curve and the parametric equation trajectory of the ideal viscoelastic elliptic geometry are read. Normalization is then performed first, and the maximum stress amplitude within the current loop is obtained. and maximum strain amplitude Divide the stress values ​​of all sampling points by Strain value divided by This generates a dimensionless point set whose values ​​are all in the range [-1, 1]. and Set the resampling arc length interval The normalized perimeter of the ideal ellipse Based on this interval, equal arc length interpolation is performed on the two dimensionless curves to construct a dimension of Given a distance matching matrix, iterate through all point pairs in the two point sets and calculate the normalized points. and The dimensionless Euclidean distance between the two curves is determined by constructing a path search tree using the discrete Fraser distance algorithm. A monotonic path that minimizes the maximum bottleneck distance is found from the lower left corner to the upper right corner of the matrix. This bottleneck distance represents the maximum topological deviation between the two curves in the dimensionless shape space. This process effectively avoids calculation errors caused by different units and ensures the objectivity of shape evaluation. Finally, the minimized maximum distance value is extracted to generate the Fraser distance difference value.

[0033] In the formula for calculating the geometric distortion index of hysteresis loops, a dimensionless index highly sensitive to shape distortion is constructed by introducing normalized topological distance and area difference. The first part... The deviation of the geometric path is quantified by the ratio of the dimensionless Fréchet distance to the perimeter. The latter part uses an exponential function to amplify the abnormal fluctuations of the energy dissipation area (hysteresis loop area), thereby enabling the identification of small nonlinear damage generated in the asphalt pavement material in the early stage of fatigue. The Frescher distance difference value represents the minimum and maximum topological distance between the closed hysteresis loop curve and the ideal viscoelastic elliptic geometry. The steps to obtain this value are as follows: Based on the dimensionless coordinate system normalized in the previous steps, calculate the optimal matching path distance between the two curves. Since it has undergone normalization, this value is dimensionless. For example, it can be calculated using an algorithm. ; To obtain the arc length of an ideal viscoelastic elliptic geometry, the steps are as follows: Read the dimensionless major axis radius of the ideal ellipse after normalization. and dimensionless minor axis radius Since the normalized ellipse is usually close to the unit circle, for example (Corresponding to the normalized value of maximum stress) (The relative width after normalization to the maximum strain reflects the phase angle characteristics.) The dimensionless perimeter is calculated using the Ramanujan approximation formula. The formula is as follows: Substituting into the calculation, we get ; The area enclosed by the closed hysteresis loop curve is obtained by: using the normalized coordinates of the closed hysteresis loop curve. The dimensionless area enclosed by the shoelace formula is calculated, which reflects the normalized energy dissipation characteristics. For example, the calculated area is... ; The area enclosed by an ideal viscoelastic ellipse is obtained by calculating the area of ​​the ideal ellipse based on the normalized major and minor axes using the following formula: Substituting the numerical values, we get .

[0034] Calculations based on parameters: Substitute the dimensionless parameters mentioned above into the main formula for calculation: First, calculate the area difference term: Index Term ; Calculate the value of the exponential function: ; Calculate the shape difference term: ; Calculate the final distortion index: ; The result indicates that the geometric distortion index of the hysteresis loop under the current loading cycle is 0.0100. This value serves as a dimensionless benchmark for subsequent assessment of the evolution trend of fatigue damage.

[0035] The steps to obtain the inflection point of the end of compressive fatigue life are as follows: Based on the hysteresis loop geometric distortion index, the hysteresis loop geometric distortion index of consecutive cycles is differentially divided according to the number of loading cycles to obtain the growth rate slope sequence of the hysteresis loop geometric distortion index. When the growth rate of the hysteresis loop geometric distortion index first exceeds the preset damage threshold, and the rate remains above the preset damage threshold in the next two consecutive cycles, the coordinate point corresponding to the cycle is determined as the inflection point of the end of the compressive fatigue life.

[0036] Specifically, based on the hysteresis loop geometric distortion index, a model is constructed that varies with the number of loading cycles. Real-time monitoring sequence of changes, setting the sliding window size Calculate the average rate of change of the distortion index within the current window as the growth rate slope. Set a preset damage threshold The threshold is set based on: accelerated loading fatigue tests on asphalt mixture specimens from the same batch, recording the distortion exponent growth rate when the specimens enter the nonlinear damage stage (accelerated point of stiffness modulus decay), and statistically calculating the average value of this rate. and standard deviation ,set up For example, as determined , ,but During real-time monitoring, the distortion exponential difference between adjacent cycles is calculated to obtain the instantaneous rate, which is then compared with... By comparison, if the current growth rate is found... Exceed And the subsequent and If all values ​​remain above the threshold, eliminating single-point mutations caused by sensor noise, then irreversible fatigue damage accumulation has occurred inside the material. The cycle number corresponding to this moment is marked as the critical point, generating the inflection point of the end of the compressive fatigue life.

Claims

1. An asphalt pavement compressive strength analysis system based on dynamic simulation, characterized in that, The system includes: The tire load discretization module is used to perform contact mechanics pre-calculation based on the three-dimensional scanning geometric data of the tire and map it to generate a ground pressure distribution cloud map. It extracts the node pressure values ​​in the ground pressure distribution cloud map, arranges the node pressure values ​​according to the geometric features of the tire tread pattern, and generates a high-resolution pressure matrix. The dynamic pressure field loading module is used to distribute pressure values ​​to the road surface finite element mesh nodes according to the high-resolution pressure matrix to generate nodal force load vectors, calculate the load transfer ratio between adjacent mesh nodes according to the nodal force load vectors, and generate a smooth dynamic load sequence. The viscoelastic hysteresis loop extraction module is used to perform time-step integral calculations on the viscoelastic constitutive parameters of asphalt pavement according to the smooth dynamic load sequence, draw and generate a closed hysteresis loop curve, and perform least squares fitting on the closed hysteresis loop curve in combination with the current dynamic modulus and phase angle parameters to generate an ideal viscoelastic elliptical geometry. The compressive fatigue life prediction module is used to calculate the Frescher distance difference value on the geometric topology based on the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry, obtain the hysteresis loop geometric distortion index, calculate the slope of the growth rate of the hysteresis loop geometric distortion index relative to the number of loading cycles, locate the coordinates of the slope change, and generate the inflection point of the end of compressive fatigue life.

2. The asphalt pavement compressive strength analysis system based on dynamic simulation according to claim 1, characterized in that, The steps for obtaining the high-resolution pressure matrix are as follows: Based on the three-dimensional scanning geometric data of the tire, the boundaries of the tread blocks and polygonal patches are analyzed, the scanning coordinate axis direction is corrected, the patch with the normal vector facing the ground is located, the contact units are aggregated according to the ground projection, the pressure values ​​are distributed according to the contact ratio of the patch, and a ground pressure distribution cloud map is generated. Based on the ground pressure distribution cloud map, the minimum and maximum values ​​of the color level mapping table are read, and the color level value of each pixel is converted into pressure units according to the corresponding index from the pixel center coordinate to the ground coordinate. Non-contact pixels are eliminated and written in the order of grid number to generate node pressure values. Based on the node pressure values, an index table is constructed according to the block numbering order of the tire tread geometry and the groove centerline sequence. First, it is sorted radially from the center of the tire crown to the tire shoulder, and then sorted circumferentially clockwise. The node pressure values ​​of each row and column are filled in sequentially to generate a high-resolution pressure matrix.

3. The asphalt pavement compressive strength analysis system based on dynamic simulation according to claim 1, characterized in that, The steps for obtaining the nodal force load vector are as follows: Based on the high-resolution pressure matrix, the vehicle's center of gravity position and roll angle parameters are read, and a projection index from the tire contact point to the road surface finite element mesh node is established. The pressure values ​​are then distributed to the road surface finite element mesh node using bilinear interpolation mapping to generate the nodal force load vector.

4. The asphalt pavement compressive strength analysis system based on dynamic simulation according to claim 1, characterized in that, The steps for obtaining the smooth dynamic load sequence are as follows: Based on the nodal force load vector, read the adjacent grid node pairs according to the grid topology, calculate the pressure difference of each pair of nodes and the ratio of the total pressure as a local proportion, write it according to the grid edge number and unify the range to generate the load transfer proportion; Based on the load transfer ratio, the nodal force load vector is weighted and distributed among adjacent grid nodes according to the time step index. The load transfer ratio is used to smooth the abrupt changes between adjacent time steps and rearrange the node order. The nodal force load vectors of each time step are then spliced ​​together to generate a smooth dynamic load sequence.

5. The asphalt pavement compressive strength analysis system based on dynamic simulation according to claim 1, characterized in that, The steps for obtaining the closed hysteresis loop curve are as follows: Based on the smooth dynamic load sequence, the time step and start and end time of the viscoelastic constitutive parameters of the asphalt pavement are read. The load increment of each time step is calculated according to the time index and the stress increment and strain increment are updated synchronously. The stress and strain values ​​are recorded in the order of the integration point number and the cycle number is marked to generate a sequence of stress and strain values ​​at key integration points. Based on the stress and strain numerical sequence of the key integration points, the start and end points of each cycle are located according to the cycle number and the corresponding intervals are extracted. The stress and strain values ​​are paired one by one according to the same time index and the curves are connected in chronological order. Discrete points outside the cycle are eliminated and the beginning and end are connected to meet the closure condition, thus generating a closed hysteresis loop curve.

6. The asphalt pavement compressive strength analysis system based on dynamic simulation according to claim 1, characterized in that, The steps for obtaining the ideal viscoelastic elliptic geometry are as follows: Based on the closed hysteresis loop curve, the dynamic modulus parameter and phase angle parameter are called to set the elliptical constraint coordinate system and initialize the major axis direction and minor axis direction. The deviation to the constraint coordinate system is calculated segment by segment according to the curve and the segment weight is set according to the dynamic modulus parameter. The major axis length, minor axis length and rotation angle are updated round by round until the total deviation no longer decreases, and an ideal viscoelastic elliptical geometry is generated.

7. The asphalt pavement compressive strength analysis system based on dynamic simulation according to claim 1, characterized in that, The steps for obtaining the hysteresis loop geometric distortion index are as follows: Based on the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry, the closed hysteresis loop curve and the ideal viscoelastic elliptical geometry are synchronously resampled at arc length intervals to establish a matching path table. The coordinates of the corresponding sampling points are matched step by step and the Euclidean distance between each pair of matching points is calculated. The minimum value of the maximum point distance is extracted along the path to generate the Fraser distance difference value. The hysteresis loop geometric distortion index is calculated based on the Fraser distance difference value.

8. The asphalt pavement compressive strength analysis system based on dynamic simulation according to claim 1, characterized in that, The steps for obtaining the inflection point of the end of the compressive fatigue life are as follows: Based on the hysteresis loop geometric distortion index, the hysteresis loop geometric distortion index of consecutive cycles is differentially divided according to the number of loading cycles to obtain the growth rate slope sequence of the hysteresis loop geometric distortion index. When the growth rate of the hysteresis loop geometric distortion index first exceeds the preset damage threshold, and the rate remains above the preset damage threshold in the next two consecutive cycles, the coordinate point corresponding to the cycle is determined as the inflection point of the end of the compressive fatigue life.

Images