A method for simulating thermal damage evolution and evaluating toughness of road material by fusing multi-scale fractal features and digital twinning

Abstract

The application provides a road material thermal damage evolution simulation and toughness evaluation method fusing multi-scale fractal features and digital twinning, comprising: collecting surface temperature change data of road material under different heating rates in real time, identifying a thermal stress concentration area in combination with a preset temperature threshold interval to obtain an initial thermal stress distribution atlas; analyzing a bimodal distance of a pore distribution curve through multi-scale fractal dimension and pore size distribution density, marking a fusion path high-risk area when the bimodal distance is lower than a preset threshold value to obtain a pore fusion trend graph; and feeding back an association database between a heating rate, multi-scale fractal dimension and damage degree through a final toughness evaluation result and a damage evolution path model to update the optimization boundary conditions of subsequent thermal damage evolution simulation.

Description

Technical Field

[0001] This invention relates to the field of information technology, and in particular to a method for simulating the thermal damage evolution and toughness assessment of road materials by integrating multi-scale fractal features and digital twins. Background Technology

[0002] Thermal damage to road materials under high-temperature environments directly affects the long-term safety and service life of pavement structures. Especially under extreme conditions such as sustained high temperatures in summer or fires, this damage evolution process significantly weakens the load-bearing capacity and durability of materials. Therefore, accurately simulating and evaluating the thermal damage toughness of road materials has become an important research direction for ensuring the reliability of transportation infrastructure. Current research generally assumes that the faster the heating rate, the greater the internal temperature gradient of the material, and the stronger the thermal stress, leading to faster pore and crack propagation, ultimately resulting in more severe damage. However, actual observations have revealed that this is not always the case. When the heating rate accelerates, small pores within the material rapidly expand and interconnect, even skipping the medium-sized pore stage and merging to form a few large pores. This change in pore size distribution from relatively uniform to highly concentrated significantly alters the damage manifestation.

[0003] This phenomenon of size transition during pore fusion presents a core technical contradiction: while an accelerated heating rate increases the overall pore volume, it also leads to an abnormal decrease in the number of medium-sized pores, while the number of large pores becomes relatively concentrated. Consequently, the number of stress concentration points decreases, and the material may retain a higher residual load-bearing capacity in the later stages of damage. This path of small pores directly fusing into large pores under rapid heating conditions contradicts the traditional understanding that damage increases monotonically with heating rate. This makes it difficult to determine which heating rate will lead to more severe material failure, and it is also impossible to reliably predict the final severity of damage based solely on the heating rate.

[0004] For example, during the heating process of asphalt mixtures, if the temperature is increased slowly at a rate of 5°C per minute, the pores gradually expand uniformly from small to medium to large, and the damage development is relatively predictable. However, when the heating rate increases to over 20°C per minute, a large number of tiny pores expand dramatically and merge rapidly in a short period of time, bypassing the medium-sized pore stage and directly forming several significantly larger pores. At this point, although the total pore area increases, the crack propagation path is restricted because the number of large pores is small and their distribution is concentrated. The material as a whole does not lose its load-bearing capacity as expected, but instead exhibits higher residual strength under certain loading conditions. This drastic and uneven change in pore size distribution caused by the heating rate renders traditional damage assessment indicators inaccurate.

[0005] Therefore, how to accurately capture the specific path of pore fusion from small to large size while considering the influence of heating rate, and how to reasonably judge the true severity of thermal damage and the toughness level of materials based on this dynamic non-uniformity of pore size distribution, has become a key issue in the simulation of thermal damage evolution and toughness assessment of road materials. Summary of the Invention

[0006] This invention provides a method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins, including:

[0007] Real-time data on surface temperature changes of road materials at different heating rates are collected, and thermal stress concentration areas are identified by combining the data with preset temperature threshold ranges to obtain an initial thermal stress distribution map.

[0008] Based on the initial thermal stress distribution map, a CT scan was used to perform a three-dimensional scan of the material interior to extract the spatial distribution density and position coordinates of small, medium and large pores, and to determine the multi-scale fractal dimension and pore size distribution density.

[0009] The spacing between the two peaks of the pore distribution curve is analyzed by the multi-scale fractal dimension and the pore size distribution density. When the spacing between the two peaks is lower than a preset threshold, it is marked as a high-risk area of ​​the fusion path, and a pore fusion trend map is obtained.

[0010] By using the pore fusion trend map, the abnormal decrease in the number of pores in the middle pore interval is extracted, and the situation where small pores directly fuse across the middle pores to form large pores is determined. Based on this, the damage evaluation benchmark is dynamically adjusted, and the distribution density of stress concentration points is determined.

[0011] Based on the stress concentration point distribution density, the load-bearing capacity attenuation gradient of the material at different temperature stages is extracted, and the multi-scale fractal influence coefficient of thermal damage on material toughness is determined by the correlation between the attenuation gradient and the multi-scale fractal dimension.

[0012] Using the multi-scale fractal influence coefficient as a damage weighting factor, a damage evolution path model is constructed in the digital twin platform. When the stress concentration point density in the high-risk area of ​​the fusion path continues to increase, the simulation boundary conditions are adjusted to obtain the toughness assessment results.

[0013] By using the toughness assessment results and the damage evolution path model, the correlation database between the heating rate, the multi-scale fractal dimension, and the degree of damage is updated to determine the optimal boundary conditions for subsequent thermal damage evolution simulation.

[0014] Furthermore, obtaining the initial thermal stress distribution map includes:

[0015] Infrared thermal imager was used to collect the surface temperature of road material specimens frame by frame. The surface temperature value of each frame was recorded according to the heating rate setting to generate a surface temperature time series curve.

[0016] The temporal temperature difference between adjacent frames is calculated based on the surface temperature time-series curve to obtain a temporal temperature difference sequence.

[0017] Extract continuous frame segments from the time-series temperature difference sequence where the time-series temperature difference value continuously exceeds a preset frame number threshold and is higher than the temperature threshold range, and mark the surface spatial location corresponding to the continuous frame segments as candidate regions for thermal stress concentration.

[0018] For the candidate region of thermal stress concentration, temperature values ​​of sampling points are extracted along the horizontal and vertical directions, and the differences at each point are calculated to form gradient distribution characteristics;

[0019] When the difference between horizontal and vertical points in the gradient distribution features both exceed the temperature threshold range, the corresponding candidate region is marked as a high thermal stress concentration region.

[0020] Based on the spatial coordinates and gradient distribution characteristics of the high thermal stress concentration area, the average value of the point-by-point difference is taken on the two-dimensional coordinate plane, and the thermal stress intensity level is divided according to the average value and mapped to the corresponding position to obtain the initial thermal stress distribution map.

[0021] Furthermore, determining the multi-scale fractal dimension and aperture distribution density includes:

[0022] Based on the spatial coordinates of the high thermal stress concentration area in the initial thermal stress distribution map, a layer-by-layer three-dimensional scan was performed using an industrial CT scanner to obtain voxel slices.

[0023] Voxels with gray values ​​below the pore gray value threshold are marked, and their three-dimensional spatial coordinates and voxel volumes are extracted to form the original pore voxel set.

[0024] For the original set of pore voxels, the equivalent diameter of the pores is calculated based on the number of voxels, and the pores are divided into small pore interval, medium pore interval and large pore interval according to the pore diameter threshold.

[0025] Using voxel grid cells as statistical windows, the proportion of pores at each scale level is statistically analyzed window by window to obtain the pore window density ratio.

[0026] The box counting method was used to process the pores at each scale level, the number of non-empty grids was counted step by step, and linear regression was performed. The absolute value of the slope of the fitted line was used as the fractal dimension of the corresponding scale, and the combination was used to obtain the multi-scale fractal dimension.

[0027] The pore window density ratios at each scale level are aligned and summarized according to three-dimensional spatial coordinates to form the pore size distribution density.

[0028] Furthermore, obtaining the pore fusion trend map includes:

[0029] Based on the pore size distribution density, an overall pore distribution curve is plotted with the equivalent pore diameter as the horizontal axis and the pore frequency of the corresponding diameter interval as the vertical axis.

[0030] Locate the local maxima of frequency on the pore distribution curve, mark the local maxima with smaller horizontal axis as small pore peaks, and mark the local maxima with larger horizontal axis as large pore peaks to obtain the bi-peak positioning coordinates;

[0031] Calculate the difference in the equivalent pore diameter corresponding to the peak value of the small pore and the peak value of the large pore to obtain the original value of the double peak spacing;

[0032] Extract the fractal dimension of the scale level to which the small hole peak belongs and the fractal dimension of the scale level to which the large hole peak belongs from the multi-scale fractal dimension, and divide the fractal dimension of the large hole peak to the fractal dimension of the small hole peak to obtain the fractal ratio value.

[0033] Multiplying the fractal ratio by the original bimodal spacing value yields the fractal-corrected bimodal spacing.

[0034] When the bimodal spacing after fractal correction is lower than the bimodal spacing threshold, the corresponding three-dimensional spatial location is marked as a high-risk area of ​​the fusion path;

[0035] Threshold comparison and marking are performed on all spatial regions covered by the pore distribution curve one by one, and all high-risk areas of the fusion path are mapped according to three-dimensional spatial coordinates to obtain a pore fusion trend map.

[0036] Furthermore, determining the stress concentration point distribution density includes:

[0037] By using the spatial location of the high-risk area of ​​the fusion path in the pore fusion trend map, the current pore frequency of the middle pore interval in that area and the pore frequency of the middle pore interval at the previous acquisition time are extracted from the pore size distribution density, and the change in the middle pore frequency is calculated.

[0038] Simultaneously, the current pore frequency of the large-pore interval in the high-risk area is extracted and compared with the pore frequency of the large-pore interval at the previous acquisition time, and the large-pore frequency increment is calculated.

[0039] When the change in the frequency of the central hole is negative and the absolute value exceeds the threshold of the sudden drop, the corresponding area is marked as the abnormal decrease area of ​​the central hole.

[0040] For the abnormal decrease region of the middle hole, the frequency increment of the large hole is compared with the absolute value of the frequency change of the middle hole. When the frequency increment of the large hole reaches a preset ratio threshold, it is determined that the small hole directly merges across the middle hole to form a large hole, and the region is marked as a cross-level fusion region.

[0041] Within the cross-level fusion region, the damage evaluation benchmark is adjusted to the cross-level fusion extension of pores;

[0042] Based on the spatial range of the cross-level fusion region, the three-dimensional spatial coordinates of each pore in the large-hole interval are located, and the coordinate position of each large hole is marked as a stress concentration point.

[0043] The total number of stress concentration points at each acquisition time within the high-risk area of ​​the fusion path is counted, and then divided by the volume of the spatial range to determine the stress concentration point distribution density sequence.

[0044] Furthermore, the multi-scale fractal influence coefficient of thermally induced damage on material toughness is determined by the correlation between the attenuation gradient and the multi-scale fractal dimension, including:

[0045] Based on the stress concentration point distribution density, loads were applied to the material specimens at each acquisition time, and the bearing capacity response values ​​were recorded.

[0046] The difference in bearing capacity response values ​​between adjacent acquisition times is calculated and divided by the corresponding temperature difference to obtain the bearing capacity attenuation gradient, thus forming an attenuation gradient sequence.

[0047] Extract the acquisition time of the abrupt change in the attenuation gradient value from the attenuation gradient sequence, and pair the stress concentration point distribution density at that time with the attenuation gradient value;

[0048] Linear regression was performed with stress concentration point distribution density as independent variable and attenuation gradient value as dependent variable, and the slope of the fitted line was used as the attenuation correlation slope.

[0049] The fractal dimensions of the small pores, medium pores, and large pores at the end of the heating process are extracted from the multi-scale fractal dimension and multiplied by the attenuation correlation slope to obtain the degradation contribution values ​​of the small pores, medium pores, and large pores.

[0050] The overall load-bearing capacity degradation level is obtained by adding the three degradation contribution values;

[0051] The overall load-bearing capacity degradation degree is divided by the initial load-bearing capacity response value before heating to obtain the load-bearing capacity degradation ratio. The load-bearing capacity degradation ratio is multiplied by the mean of the small-pore fractal dimension, medium-pore fractal dimension and large-pore fractal dimension to determine the multi-scale fractal influence coefficient.

[0052] Furthermore, the obtained toughness assessment results include:

[0053] The multi-scale fractal influence coefficient is used as the damage weighting factor to establish a three-dimensional finite element mesh in the digital twin platform;

[0054] Map the three-dimensional spatial coordinates of the high-risk areas of the fusion path in the porosity fusion trend diagram to the corresponding elements of the finite element mesh;

[0055] For the mapped mesh element, the preset elastic modulus is divided by the damage weighting factor to obtain the reduced elastic modulus, thus forming a damage-weighted finite element mesh.

[0056] Initial simulation boundary conditions are set, including temperature load and mechanical load. The initial simulation boundary conditions are used to drive the finite element mesh to solve for the temperature field and stress field, and the initial simulation results of damage evolution are obtained.

[0057] At each simulation moment, the density of simulated stress concentration points in the high-risk area of ​​the fusion path is statistically analyzed. When the density of simulated stress concentration points continues to increase for multiple consecutive simulation moments and the increase exceeds the density growth threshold, the simulation boundary conditions are adjusted, additional temperature load increments are added, and the finite element mesh is re-driven to obtain updated damage evolution simulation results.

[0058] Based on the updated damage evolution simulation results, the residual bearing capacity values ​​of all grid elements at the final simulation time are extracted, the residual bearing capacity ratio is calculated, and the residual bearing capacity ratio is divided by the damage weight factor to obtain the toughness assessment results.

[0059] Furthermore, determining the optimized boundary conditions for subsequent thermally induced damage evolution simulation includes:

[0060] The current heating rate level, the multi-scale fractal dimension, the overall load-bearing capacity degradation degree, and the toughness assessment result are recorded as a single associated record and written into the associated database.

[0061] When a record with the same heating rate already exists in the associated database, the existing record is replaced with the current associated record to complete the feedback update;

[0062] Based on the updated associated database, retrieve the overall load-bearing capacity degradation and toughness assessment results corresponding to all heating rate levels;

[0063] Match the heating rate setting with the highest overall load-bearing capacity degradation value, extract the temperature load and load adjustment coefficient used for the heating rate setting, and determine the optimized boundary conditions for subsequent thermal damage evolution simulation.

[0064] The technical solutions provided by the embodiments of the present invention may include the following beneficial effects:

[0065] This invention discloses a method for simulating the thermal damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins. The method acquires real-time surface temperature change data of road materials at different heating rates, identifies thermal stress concentration areas by combining preset temperature threshold ranges, and generates an initial thermal stress distribution map. Then, CT scans are used to extract the spatial distribution density and location coordinates of small, medium, and large pores within the material, calculating the multi-scale fractal dimension and pore size distribution density. High-risk areas of the fusion path are marked when the bimodal distance of the pore distribution curve is below a threshold, and the damage evaluation benchmark is dynamically adjusted based on the anomalous drop in the number of medium pores to accurately locate the stress concentration point distribution density. Furthermore, the bearing capacity attenuation gradient at different temperature stages is extracted and correlated with the multi-scale fractal dimension to quantify the multi-scale fractal influence coefficient of thermal damage on material toughness. Using this coefficient as a damage weighting factor, a damage evolution path model is constructed in the digital twin platform, dynamically adjusting the simulation boundary conditions to achieve toughness assessment. Finally, the assessment results and the evolution model are fed back to update the database of heating rate, fractal dimension, and damage degree correlation, optimizing the subsequent thermal damage simulation boundary conditions. This invention solves the core problem of rapid toughness degradation in road materials caused by thermal stress-induced multi-scale pore fusion under complex thermal environments, significantly improving the accuracy of thermal damage evolution prediction and the reliability of material toughness assessment, and providing scientific support for the durability design and maintenance of road materials. Attached Figure Description

[0066] Figure 1 This is a flowchart of a method for simulating the thermal damage evolution and toughness assessment of road materials that integrates multi-scale fractal features and digital twins, according to the present invention.

[0067] Figure 2 This is a schematic diagram of a method for simulating the thermal damage evolution and assessing the toughness of road materials that integrates multi-scale fractal features and digital twins, according to the present invention.

[0068] Figure 3 This is another schematic diagram of a method for simulating the thermal damage evolution and toughness assessment of road materials that integrates multi-scale fractal features and digital twins according to the present invention.

[0069] Figure 4 This is a display of the equivalent diameter of each temperature orifice in an embodiment of the present invention.

[0070] Figure 5 The internal pore structure at various temperatures in the embodiments of the present invention is shown in (a)-(d) as a schematic diagram of the reconstruction of temperature models at 25°C, 250°C, 500°C, and 750°C. (e) a comparison of the pore structure between the normal temperature model and the 1000°C model. Detailed Implementation

[0071] To further understand the content of this invention, a detailed description of the invention is provided in conjunction with the accompanying drawings and embodiments. The specific embodiments described herein are for illustrative purposes only and are not intended to limit the invention. It should also be noted that, for ease of description, only the parts relevant to the invention are shown in the accompanying drawings.

[0072] like Figure 1 As shown in the figure, this embodiment of a method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins may specifically include:

[0073] Step S101: Real-time acquisition of surface temperature change data of road materials at different heating rates, combined with preset temperature threshold range to identify areas of concentrated thermal stress, and obtain initial thermal stress distribution map.

[0074] Infrared thermal imagers are used to collect frame-by-frame temperature data on the surface of road material specimens. The surface temperature value for each frame is recorded according to a preset heating rate level. The collected frame-by-frame surface temperature values ​​are arranged chronologically to generate a surface temperature time-series curve. Based on the temperature value changes between adjacent frames in the surface temperature time-series curve, the time-series temperature difference is calculated frame-by-frame to obtain a time-series temperature difference sequence corresponding to the heating rate level. Based on the time-series temperature difference sequence, consecutive frame segments whose time-series temperature difference values ​​continuously exceed a preset frame number threshold and remain above a preset temperature threshold range are extracted. The surface spatial locations corresponding to these consecutive frame segments are marked as candidate regions for thermal stress concentration. For each candidate region, temperature values ​​at each sampling point are extracted along the horizontal and vertical directions of the material surface. The point-by-point difference between the horizontal and vertical sampling point temperatures is calculated. The spatial arrangement of these point-by-point differences at each sampling point is used as the gradient distribution feature of the candidate region. If both the horizontal and vertical point-by-point differences in the gradient distribution feature simultaneously exceed a preset temperature threshold range, the candidate region is marked as a high thermal stress concentration region. Based on the spatial coordinates and corresponding gradient distribution characteristics of the high thermal stress concentration area, the average difference of each sampling point in each high thermal stress concentration area is taken on the two-dimensional coordinate plane of the material surface. The thermal stress intensity level is divided according to the magnitude of the average value, and the thermal stress intensity level is mapped to the corresponding position on the two-dimensional coordinate plane in the form of regional markers to obtain the initial thermal stress distribution map.

[0075] In this embodiment, the temperature threshold range for asphalt mixtures is set to 60-80℃, corresponding to the temperature range near the softening point of asphalt; for cement mortar or concrete, the temperature threshold range is set to 150-300℃, corresponding to the temperature range of internal moisture evaporation and Ca(OH)2 decomposition in cement-based materials. This threshold is typically determined through material thermal analysis experiments, such as differential scanning calorimetry (DSC).

[0076] In one embodiment, the infrared thermal imager continuously scans the entire exposed surface of the road material specimen at a fixed frame rate. The surface temperature value collected in each frame covers all pixel positions on the specimen surface. The heating rate settings are several sets of constant heating rates preset on the heating device. Each set of heating rate settings corresponds to an independent frame-by-frame acquisition process. The surface temperature values ​​of each frame are arranged in chronological order of acquisition time to form the surface temperature time sequence curve.

[0077] Specifically, the time-series temperature difference refers to the temperature difference between two adjacent frames at the same pixel location in the surface temperature time-series curve. After traversing all pixels frame by frame, each frame corresponds to a set of time-series temperature difference values. The time-series temperature difference values ​​of all frames are summarized to form the time-series temperature difference sequence.

[0078] In one possible implementation, the criteria for determining the consecutive frame segments include two conditions: first, the temporal temperature difference value of each frame in the segment is higher than the preset temperature threshold range; second, the number of consecutive frames that meet the temperature difference condition reaches the preset frame number threshold.

[0079] In this embodiment of the application, the preset frame rate threshold is determined according to the acquisition frame rate of the infrared thermal imager, and is usually set to 5-10 consecutive frames. When the acquisition frame rate is 25fps, it corresponds to a duration of 0.2-0.4 seconds, ensuring that the detected temperature anomaly is a continuous change rather than instantaneous noise.

[0080] For example, if the temporal temperature difference of a certain pixel location exceeds the preset temperature threshold range for several consecutive frames and the number of frames reaches the preset frame number threshold, then the surface spatial location of the pixel is marked as a candidate region for high thermal stress concentration. This is the initial screening stage. After further screening, it can be confirmed as a region for high thermal stress concentration. This calibration process is performed on all pixels on the surface of the specimen one by one.

[0081] It should be noted that the gradient distribution characteristics are obtained as follows: Within each candidate region of thermal stress concentration, adjacent sampling points are selected sequentially along the horizontal direction, and the difference in temperature values ​​between two adjacent sampling points is calculated. After traversing all pairs of adjacent sampling points in the horizontal direction, a horizontal point-by-point difference sequence for the region is obtained. The same operation is performed along the vertical direction to obtain a vertical point-by-point difference sequence. Each difference in the horizontal and vertical point-by-point difference sequences corresponds to the spatial coordinates of a sampling point. Arranging all differences according to the spatial coordinates of the corresponding sampling points constitutes the gradient distribution characteristics of the candidate region. The temperature threshold interval is a preset temperature gradient critical value based on the material type. If the horizontal point-by-point difference and the vertical point-by-point difference in the gradient distribution characteristics of a candidate region simultaneously exceed the temperature threshold interval at the same sampling point, then the candidate region is marked as a high thermal stress concentration region.

[0082] In one embodiment, the thermal stress intensity level is divided as follows: the arithmetic mean of the point-by-point differences of all sampling points in each high thermal stress concentration area is taken, and a corresponding level label is assigned according to the preset numerical range in which the mean falls, wherein the numerical range corresponds to the thermal stress intensity from low to high level in order from low to high.

[0083] Preferably, the initial thermal stress distribution map uses a two-dimensional coordinate plane of the material surface as the base map. Each high thermal stress concentration area is marked with a region marker indicating its thermal stress intensity level at the corresponding position on the base map. Unmarked areas indicate that the thermal stress concentration judgment has not been triggered at the current heating rate setting, thereby forming a complete initial thermal stress distribution map.

[0084] Step S102: Based on the initial thermal stress distribution map, a CT scan is used to perform a three-dimensional scan of the material interior, extract the spatial distribution density and position coordinates of small, medium and large pores, and determine the multi-scale fractal dimension and pore size distribution density.

[0085] Based on the spatial coordinates of the high thermal stress concentration area in the initial thermal stress distribution map, an industrial CT scanner is used to perform layer-by-layer three-dimensional scanning of the road material specimen to obtain all voxel slices covering the high thermal stress concentration area. Voxels with gray values ​​lower than a preset pore gray value threshold are marked, and the three-dimensional spatial coordinates and voxel volumes of the marked voxels are extracted to obtain the original pore voxel set. For each pore in the original pore voxel set, the equivalent pore diameter is calculated based on the number of voxels it occupies. According to a preset pore diameter threshold, the equivalent pore diameter is divided into three scale levels: small pore, medium pore, and large pore. For each scale level, a preset voxel grid unit is used as a statistical window in three-dimensional space to count the ratio of the number of pores of that scale level to the total number of voxels in the window, obtaining the pore window density ratio for each scale level. Based on the pore space distribution corresponding to the pore window density ratios at each scale level, box counting is used to process the small-pore, medium-pore, and large-pore intervals respectively: A progressively smaller cubic mesh is used to cover the three-dimensional space containing the pores at each scale level. The number of non-empty meshes containing at least one pore voxel is counted at each level. The logarithm of the mesh side length and the corresponding number of non-empty meshes at each level is then used for linear regression fitting. The absolute value of the slope of the fitted line is the fractal dimension at that scale level. The fractal dimensions of small-pore, medium-pore, and large-pore are obtained respectively, and combined to determine the multi-scale fractal dimension. Using the pore window density ratios at each scale level, the window density ratios of the small-pore, medium-pore, and large-pore intervals are aligned and summarized according to three-dimensional spatial coordinates to form a spatial density distribution covering all scale levels, thus determining the pore size distribution density.

[0086] In this embodiment, the pore size threshold is calibrated according to the type of road material. Specifically, for asphalt mixtures, based on the interface characteristics between asphalt mastic and aggregates, the small pore range is defined as having an equivalent pore diameter < 50 μm (corresponding to asphalt film thickness and microbubbles), the medium pore range as 50-200 μm (corresponding to pores at the aggregate-asphalt interface), and the large pore range as ≥ 200 μm (corresponding to interconnected pores between aggregates). For cement mortar / concrete: the small pore range is defined as having an equivalent pore diameter < 100 μm (corresponding to gel pores and capillary pores), the medium pore range as 100-500 μm (corresponding to pores in the aggregate interface transition zone), and the large pore range as ≥ 500 μm (corresponding to pores and construction defect pores). The above thresholds are determined based on the frequency valley position of the initial pore size distribution curve of the material, the microstructural characteristic scale of the material, and the classification criteria for pores in industry standards. The threshold boundaries can be adaptively determined based on the CT scan results of the specific material and the valley position of the pore size distribution curve.

[0087] In one embodiment, the industrial CT scanner scans the road material specimen layer by layer according to a preset scanning layer spacing. Each layer scan produces a voxel slice, and each voxel in each voxel slice records the gray value of that location and its three-dimensional spatial coordinates inside the specimen. The scanning range covers the specimen depth range corresponding to all high thermal stress concentration areas in the initial thermal stress distribution map.

[0088] Specifically, the preset porosity grayscale threshold is determined as follows: A known non-porous dense region in the specimen is taken as a reference. The average grayscale value of all voxels within this dense region is calculated. This average value is multiplied by a preset attenuation ratio to obtain the porosity grayscale threshold. In all voxel slices, voxels with grayscale values ​​below this threshold are identified as porosity voxels and marked. The three-dimensional spatial coordinates and voxel volumes of all marked voxels constitute the original porosity voxel set.

[0089] It should be noted that the conversion process for the equivalent pore diameter is performed on each individual pore in the original pore voxel set. An individual pore consists of several spatially adjacent marked voxels. Adjacent marked voxels are grouped into the same pore using a 26-neighborhood connectivity determination. The sum of the voxel volumes of all marked voxels contained in each individual pore is the actual volume of that pore. This actual volume is considered as the volume of a sphere of equal volume, and the diameter of this sphere is calculated to obtain the equivalent pore diameter. This equivalent method is a commonly used approximation in the evaluation of porosity in road materials. Although pore morphologies vary greatly, it provides a reasonable size representation in subsequent fractal dimension calculations, and potential inaccuracies can be mitigated through statistical analysis.

[0090] For example, if the total volume of the marked voxels contained in a certain independent pore corresponds to a sphere, the diameter of the sphere is the pore equivalent diameter of the pore.

[0091] In one possible implementation, the pore size threshold includes two preset diameter limit values, dividing the equivalent diameter of all pores into three scale levels: pores with equivalent diameters smaller than the first diameter limit value are classified into the small pore interval, those between the first and second diameter limit values ​​are classified into the medium pore interval, and those larger than the second diameter limit value are classified into the large pore interval. Further, the statistical process for the pore window density ratio at each scale level is as follows: In three-dimensional space, a cube with a preset side length is used as the statistical window, dividing the three-dimensional scanning space of the specimen into several non-overlapping statistical windows. For the small pore interval, each statistical window is traversed, and the number of pores belonging to the small pore interval within that window is counted. This number is then divided by the total number of elements contained within that window to obtain the small pore density ratio for that window. After performing this operation on all statistical windows, the distribution of the pore window density ratio for the small pore interval in three-dimensional space is formed. The medium pore interval and the large pore interval are statistically analyzed window by window in the same way, each forming a spatial distribution of the pore window density ratio for its corresponding scale level.

[0092] For example, the specific operation of the box counting method is performed independently for each scale level. Taking the small hole region as an example, within the three-dimensional space occupied by the pores in the small hole region, a cube mesh with an initial side length is used to fully cover the space, and the number of non-empty meshes containing at least one small hole voxel is counted and recorded. Then, the side length of the cube mesh is reduced to half of the previous level, and the same space is fully covered again, and the number of non-empty meshes is counted. In this way, the mesh side length is reduced level by level, and the corresponding mesh side length and the number of non-empty meshes are recorded at each level, until the mesh side length is reduced to the size of a single voxel.

[0093] Understandably, the grid edge lengths and the number of non-empty grids recorded at each level constitute a set of discrete data points. The natural logarithm of the grid edge length for each data point is used as the horizontal axis value, and the natural logarithm of the corresponding number of non-empty grids is used as the vertical axis value. A linear regression is then performed on all data points on a coordinate plane formed by the horizontal and vertical axis values. The absolute value of the slope of the fitted line is the fractal dimension of the pore interval. This fractal dimension reflects the complexity of the pore filling in three-dimensional space; a larger value indicates a more complex and irregular spatial distribution of pores.

[0094] Preferably, the medium-sized aperture interval and the large-sized aperture interval are executed independently according to the same process of the box counting method described above, and the fractal dimension of the medium-sized aperture and the fractal dimension of the large-sized aperture are obtained respectively. The fractal dimensions of the small-sized aperture, the medium-sized aperture, and the large-sized aperture are arranged and combined according to the scale level to determine the multi-scale fractal dimension.

[0095] In one embodiment, the pore size distribution density is determined as follows: First, the pore window density ratios for the small pore interval, medium pore interval, and large pore interval are calculated, which is the ratio of the pore volume to the statistical window volume in each interval; then, these density ratios are aligned window by window according to the three-dimensional spatial coordinates of the statistical window, and the density ratios of the three scale levels are summarized at the same statistical window position to form the comprehensive density value of the window at all scale levels. The distribution of the comprehensive density values ​​of all statistical windows in three-dimensional space constitutes the pore size distribution density.

[0096] Step S103: Analyze the bimodal spacing of the pore distribution curve using multi-scale fractal dimension and pore size distribution density. When the bimodal spacing is lower than a preset threshold, it is marked as a high-risk area for the fusion path, thus obtaining a pore fusion trend map.

[0097] like Figure 2 As shown, based on the pore size distribution density, pores at all scale levels within the statistical window are summarized according to their equivalent diameter. A pore distribution curve is plotted with the equivalent diameter as the horizontal axis and the pore frequency within the corresponding diameter range as the vertical axis. The locations of local maxima of frequency are identified on this curve. Local maxima with smaller equivalent diameters on the horizontal axis are marked as small-pore peaks, and local maxima with larger equivalent diameters on the horizontal axis are marked as large-pore peaks, thus obtaining the bimodal positioning coordinates. Based on these bimodal positioning coordinates, the difference between the equivalent diameters of the small-pore peak and the large-pore peak is calculated to obtain the original bimodal spacing. The fractal dimensions of the scale levels to which the small-pore peak belongs and the large-pore peak belongs are extracted from the multi-scale fractal dimension. The fractal ratio is obtained by dividing the fractal dimension of the large-pore peak by the fractal dimension of the small-pore peak. This fractal ratio is then multiplied by the original bimodal spacing value to obtain the fractally corrected bimodal spacing. If the bimodal spacing after fractal correction is lower than the preset bimodal spacing threshold, the three-dimensional spatial location corresponding to that region is marked as a high-risk area for fusion path. The threshold comparison and marking are performed one by one on all spatial regions covered by the pore distribution curve. All marked high-risk areas for fusion path are mapped onto the spatial distribution map inside the material according to the three-dimensional spatial coordinates to obtain a pore fusion trend map.

[0098] Specifically, in the embodiments of this application, the peak spacing threshold is determined based on the initial pore size distribution characteristics of the material, and is usually set to 30%-50% of the initial peak spacing. For example, if the peak value of the small pore corresponds to 50μm and the peak value of the large pore corresponds to 400μm in the initial state of the material, then the peak spacing threshold is set to (400-50)×0.4=140μm.

[0099] In one embodiment, the overall pore distribution curve is plotted based on the pore size distribution density of all statistical windows. Specifically, all identified pores inside the material are classified into several preset diameter intervals according to their equivalent diameter. The number of pores in each diameter interval is counted as the frequency value of that interval. Points are plotted with the diameter interval as the horizontal axis and the frequency value as the vertical axis, and a continuous curve reflecting the distribution of all pores in different diameter ranges is formed.

[0100] Specifically, the process of obtaining the bimodal positioning coordinates involves traversing the frequency values ​​point by point on the overall pore distribution curve. When the frequency value of a certain diameter interval is simultaneously greater than the frequency values ​​of its left and right adjacent intervals, that point is determined to be a local maximum. Among all local maxima, they are arranged in ascending order according to their corresponding equivalent pore diameters on the horizontal axis. Local maxima with smaller diameters are marked as small-pore peak values, and local maxima with larger diameters are marked as large-pore peak values. Their positions on the horizontal axis constitute the bimodal positioning coordinates.

[0101] It should be noted that the original value of the double peak spacing is the difference between the equivalent pore diameter value corresponding to the large pore peak value and the equivalent pore diameter value corresponding to the small pore peak value. This difference represents the absolute distance between the two peaks in the diameter dimension.

[0102] In one possible implementation, the fractal ratio is constructed and corrected as follows: The scale level of the small aperture peak is determined based on its equivalent pore diameter; the fractal dimension of that scale level is extracted from the multi-scale fractal dimension. Similarly, the scale level of the large aperture peak is determined based on its equivalent pore diameter, and its corresponding fractal dimension is extracted. The fractal dimension of the scale level of the large aperture peak is divided by the fractal dimension of the scale level of the small aperture peak to obtain the fractal ratio. This fractal ratio reflects the relative difference in spatial filling complexity between the large and small apertures. When the spatial distribution complexity of the large aperture is higher than that of the small aperture, the fractal ratio is greater than 1. In this case, the corrected bimodal spacing is amplified, indicating that the pore size span between the two peaks is further widened in the spatial structure. When the fractal ratio is less than 1, the corrected bimodal spacing is reduced, indicating that the two peaks tend to be closer in the spatial structure. Multiplying the fractal ratio by the original bimodal spacing value yields the fractally corrected bimodal spacing. Furthermore, the preset bi-peak spacing threshold is a critical spacing value pre-calibrated according to the road material type. When the bi-peak spacing after fractal correction is lower than this threshold, it indicates that the distance between the peak value of the small hole and the peak value of the large hole in the diameter dimension has been shortened to the critical range of cross-level fusion of pores. At this time, the corresponding spatial area is marked as a high-risk area of ​​fusion path.

[0103] For example, the porosity fusion trend map is formed as follows: For each spatial region within the material, the bimodal spacing after fractal correction is compared with a preset threshold. The spatial locations of all regions marked as high-risk areas for fusion paths are then mapped onto a spatial distribution map within the material using three-dimensional spatial coordinates, forming a complete porosity fusion trend map. Specifically, the formula for calculating the bimodal spacing after fractal correction is as follows: D corrected =(D large / D small )×Δd, where D corrected D is the bimodal spacing after fractal correction. large D is the fractal dimension of the scale level to which the large aperture peak belongs. small Δd is the fractal dimension of the scale level to which the pinhole peak belongs, and Δd is the original value of the bipeak spacing.

[0104] Step S104: Extract the abnormal decrease in the number of pores in the middle pore interval through the pore fusion trend map. Since the sudden drop in the number of middle pores indicates that small pores directly fuse across the middle pores to form large pores, the damage evaluation benchmark is dynamically adjusted accordingly to determine the distribution density of stress concentration points.

[0105] Based on the spatial location of the high-risk area of ​​the fusion path marked in the pore fusion trend map, the current pore frequency of the medium-sized pores in this area and the pore frequency of the medium-sized pores at the previous acquisition time in the same heating process are extracted from the pore size distribution density. The change in medium-sized pore frequency is obtained by subtracting the medium-sized pore frequency at the previous acquisition time from the current medium-sized pore frequency. At the same time, the current pore frequency of the large-sized pores in this area and the pore frequency of the large-sized pores at the previous acquisition time are extracted. The increase in large-sized pore frequency is obtained by subtracting the large-sized pore frequency at the previous acquisition time from the current large-sized pore frequency. If the change in medium-sized pore frequency is negative and its absolute value exceeds the preset threshold for sudden drop, then this area is marked as an abnormal decrease area of ​​medium-sized pores. For the abnormal decrease region of the medium pore, the frequency increment of the large pore is compared with the absolute value of the frequency change of the medium pore. If the frequency increment of the large pore reaches or exceeds the preset proportional threshold of the absolute value of the frequency change of the medium pore, it is determined that the small pores in the region directly merge across the medium pore to form a large pore. This region is marked as a cross-level fusion region. In the cross-level fusion region, the damage evaluation benchmark is adjusted from the default uniform progressive expansion of pores to cross-level fusion expansion of pores.

[0106] Specifically, the method for determining the preset ratio threshold is based on the principle of pore volume conservation. When the pores in the mesopore interval are transformed into the pores in the macropore interval, the factors affecting the preset ratio threshold include the pore closure loss rate and the interference of new pore generation.

[0107] In a more detailed embodiment, some mesopores will close due to material expansion at high temperatures, with a typical closure rate of approximately 10%-20%. New micropores may also be generated within the material at high temperatures, contributing to the increase in macropore frequency, with a typical closure rate of approximately 5%-15%. Considering these factors, the recommended range for the preset proportion threshold is 0.6-0.8. The specific value can be determined through the following calibration experiment: Perform multiple repeated CT scans on the same material under standard heating conditions, statistically analyze the correspondence between the disappearance of mesopores and the increase in macropores, and use the minimum conversion rate at a 95% confidence interval as the proportion threshold. For example, for dense-graded asphalt mixtures, the preset proportion threshold is 0.7; for cement mortar, the preset proportion threshold is 0.65.

[0108] The threshold for the magnitude of the sudden drop is used to determine whether there is an abnormal decrease in the number of mesopores. This threshold is set to 20%-30% of the mesopore frequency at the previous sampling time. For example, if the mesopore frequency at the previous time was 100, the threshold for the magnitude of the sudden drop is set to 20-30. When the change in mesopore frequency is negative and the absolute value exceeds 20-30, it is determined to be an abnormal decrease in the number of mesopores. This threshold can also be calibrated through pre-experiments based on the porosity evolution characteristics of the material.

[0109] Based on the spatial range of the cross-level fusion region, the three-dimensional spatial coordinates of each pore in the large pore interval are located one by one within the region. The coordinate position of each large pore is marked as a stress concentration point. At each acquisition time during the same heating process, the total number of stress concentration points in the entire spatial range of the high-risk area of ​​the fusion path is counted and divided by the volume of the spatial range to determine the stress concentration point distribution density sequence at each time.

[0110] In one implementation, the extraction of the frequency change of mesopores and the frequency increment of macropores is performed within a high-risk region of the fusion path marked by the pore fusion trend map. For this region, the pore frequency values ​​of the mesopore region at the current acquisition time are read from the pore size distribution density, and simultaneously the pore frequency values ​​of the macropore region at the current acquisition time are read. These values ​​are then subtracted from the corresponding frequency values ​​at the previous acquisition time to obtain the frequency change of mesopores and the frequency increment of macropores.

[0111] It should be noted that the "previous acquisition time" refers to the time during which the industrial CT scanner performs multiple consecutive three-dimensional scans on the specimen at preset scanning intervals during the same heating process. Each scan corresponds to one acquisition time, and the scan preceding the current acquisition time is the previous acquisition time. Each scan independently extracts the pore size distribution density; therefore, the same spatial region has its own corresponding mesopore frequency and macropore frequency at different acquisition times.

[0112] In a preferred embodiment of this application, the internal pores of the cement mortar are divided into four levels based on the equivalent diameter of the pores: small pores (<100μm), medium pores (100-500μm), large pores (500-1000μm), and ultra-large pores (≥1000μm). The three-dimensional pore models after each temperature treatment are then classified and extracted according to the pore radius based on the predefined pore levels. The extraction results are as follows: Figure 4 As shown in the figure, the analysis results of the pore model show that temperature treatment has a significant impact on the pore structure of mortar samples, and different temperatures have a clear selective effect on the development of pores at different scales.

[0113] Specifically, the determination process for the cross-level fusion region is as follows: Within the region of abnormal decrease in the number of medium-sized holes, the absolute value of the increase in the frequency of large holes is compared with that of the change in the frequency of medium-sized holes. This comparison is to verify whether the pores that have decreased in the medium-sized hole region have been transformed into pores in the large-sized hole region. If the change in the frequency of medium-sized holes is negative, it indicates that the number of medium-sized holes has decreased at the current acquisition time relative to the previous acquisition time. If, at the same time, the increase in the frequency of large holes is positive and reaches a preset proportional threshold above the absolute value of the change in the frequency of medium-sized holes, it indicates that the pores that have disappeared in the medium-sized hole region have not closed on their own, but have merged with the surrounding small holes and directly crossed the medium-sized hole scale to form large holes.

[0114] For example, in asphalt mixture specimens, when rapid heating causes the micropores to expand violently, adjacent small pores connect and bypass the process of gradually expanding into mesopores, directly forming large pores with significantly larger sizes. At this time, the number of mesopores decreases sharply while the number of large pores increases sharply. The area that meets the above comparison conditions is marked as the cross-level fusion area.

[0115] In one possible implementation, the adjustment of the damage assessment benchmark refers to the following: in regions where cross-level fusion has not occurred, pore evolution follows a uniform progressive path, gradually expanding from small pores to medium pores and then from medium pores to large pores; while in regions where cross-level fusion occurs, the pore evolution path skips the medium pore stage, and the damage assessment benchmark is correspondingly adjusted from uniform progressive pore expansion to cross-level fusion expansion of pores. Furthermore, the calibration of stress concentration points is performed within the cross-level fusion region. Since large pores in this region are formed by the merging of small pores across levels, the three-dimensional spatial coordinates of each large pore are calibrated as a stress concentration point, and the coordinates of all large pores constitute the spatial set of stress concentration points within this region.

[0116] Voxelstudio Recon software was opened, and CT scan data after temperature treatment was imported. Median filtering was used to reduce noise and improve the image signal-to-noise ratio. For hardening artifacts, edge hardening was removed using reconstruction software, and a three-dimensional pore model was reconstructed. The CT slices were scanned layer by layer, and the model was reconstructed. Differences in grayscale values ​​were used to distinguish the internal components of the material. Cracks and pores were separated from the reconstructed model, thresholded, and the internal pores were classified and labeled according to their size. The entire process is as follows: Figure 5 As shown, before the heat treatment, there are threshold voids in the edge region. After the heat treatment, it can be observed that the number and density of pores around the pre-fabricated crack gradually increase with the effect of the temperature gradient. As the temperature rises around the crack, the crack expands, the pores merge, and new microcracks appear.

[0117] For example, the stress concentration point distribution density is determined as follows: the total number of stress concentration points within the entire spatial range of the high-risk area of ​​the fusion path is counted, and then divided by the volume of that spatial range. The quotient is the stress concentration point distribution density, which reflects the density of stress concentration points per unit volume. Specifically, the formula for calculating the stress concentration point distribution density is as follows: ρ stress =N points / V region , where ρ stress N represents the stress concentration point distribution density. points V represents the total number of stress concentration points. region The volume of high-risk areas in the fusion path.

[0118] Step S105: Based on the stress concentration point distribution density, extract the load-bearing capacity attenuation gradient of the material at different temperature stages, evaluate the overall load-bearing performance degradation degree by the correlation between the attenuation gradient and the multi-scale fractal dimension, and determine the multi-scale fractal influence coefficient of thermal damage on the material toughness.

[0119] Based on the stress concentration point distribution density, at each sampling time during the same heating process, a preset load is applied to the material specimen and the corresponding bearing capacity response value is recorded. The difference between the bearing capacity response values ​​of two adjacent sampling times is divided by the temperature difference extracted from the surface temperature time series curve between the two sampling times to obtain the bearing capacity attenuation gradient between each pair of adjacent time times. The bearing capacity attenuation gradients corresponding to all sampling times are arranged in chronological order to obtain the attenuation gradient sequence. Through the attenuation gradient sequence, the sampling times where the attenuation gradient value changes abruptly are extracted. The stress concentration point distribution density corresponding to the time is paired with the attenuation gradient value. For the paired data of these abrupt time times, a linear regression is performed with the stress concentration point distribution density as the independent variable and the attenuation gradient value as the dependent variable. The slope of the fitted line is the attenuation correlation slope, which reflects the change in the bearing capacity attenuation gradient when the stress concentration point distribution density increases by one unit. Based on the attenuation correlation slope, the fractal dimensions of the small pore, medium pore, and large pore at the end of the heating process are extracted from the multi-scale fractal dimension. These three fractal dimensions are multiplied by the attenuation correlation slope to obtain the degradation contribution values ​​of the small pore, medium pore, and large pore, respectively. The sum of these three degradation contribution values ​​yields the overall load-bearing capacity degradation degree. This overall load-bearing capacity degradation degree is divided by the initial load-bearing capacity response value recorded when the material specimen was subjected to the same preset load before heating to obtain the load-bearing capacity degradation ratio. This ratio is then multiplied by the mean of the small pore, medium pore, and large pore fractal dimensions in the multi-scale fractal dimension to determine the multi-scale fractal influence coefficient of thermally induced damage on material toughness.

[0120] In one embodiment, the load-bearing capacity response value is obtained as follows: at each sampling moment during the same heating process, a preset constant load is applied to the surface of the material specimen, and the mechanical response reading generated by the specimen under this load is recorded. This reading is the load-bearing capacity response value at that sampling moment. The preset load remains consistent throughout all sampling moments to eliminate the interference of load differences on the response value.

[0121] Specifically, the calculation of the bearing capacity attenuation gradient is performed between two adjacent acquisition times. The bearing capacity response value at the later acquisition time is subtracted from the bearing capacity response value at the earlier acquisition time to obtain the bearing capacity change difference. Simultaneously, the surface temperature values ​​corresponding to the later and earlier acquisition times are extracted from the surface temperature time-series curves mentioned above, and the difference between the two is obtained as the temperature difference. The bearing capacity change difference is divided by the temperature difference, and the quotient is the bearing capacity attenuation gradient between the adjacent time pairs, calculated using the formula: G decay =(F {t+1} -F t ) / (T{t+1} -T t ), where G decay For the load-bearing capacity attenuation gradient, F t and F {t+1} The bearing capacity response values ​​at the previous and next data acquisition times, respectively, T t and T {t+1} These represent the corresponding surface temperature values. The load-bearing capacity decay gradient characterizes the change in load-bearing capacity response value for each unit increase in temperature.

[0122] It should be noted that the attenuation gradient sequence is composed of all adjacent time points with respect to their respective load-bearing capacity attenuation gradients arranged in chronological order of acquisition time. The changing trend of the attenuation gradient values ​​in the sequence reflects the rate of change of the material's load-bearing capacity during the heating process. Furthermore, the paired data is constructed as follows: in the attenuation gradient sequence, for each acquisition time, the stress concentration point distribution density value corresponding to that time is paired with the attenuation gradient value corresponding to that time to form a set of paired data. The paired data for all acquisition times are then aggregated to form a two-dimensional scatter plot data set.

[0123] In one possible implementation, the linear regression fitting process is as follows: using the stress concentration point distribution density as the independent variable and the attenuation gradient as the dependent variable, all data points in the aforementioned scattered dataset are marked on a two-dimensional coordinate plane, and a straight line is fitted to all data points using the least squares method. The core operation of the least squares method is to find a straight line that minimizes the sum of the squares of the vertical distances from all data points to that line. After fitting, the slope of the straight line is the attenuation correlation slope. When the attenuation correlation slope is negative, it indicates that the greater the stress concentration point distribution density, the faster the rate of decrease of the bearing capacity attenuation gradient; when the absolute value of the attenuation correlation slope is larger, it indicates that the stress concentration point distribution density has a more significant impact on bearing capacity attenuation.

[0124] For example, during the rapid heating process of asphalt mixture, if the concentration of large pores in the cross-grade fusion region leads to a rapid increase in the distribution density of stress concentration points, and the decrease in the bearing capacity attenuation gradient also intensifies, the absolute value of the attenuation correlation slope obtained by fitting will be large.

[0125] For example, the degradation contribution value is calculated as follows: the fractal dimension of the small pore is multiplied by the attenuation correlation slope to obtain the degradation contribution value of the small pore; the fractal dimension of the medium pore is multiplied by the attenuation correlation slope to obtain the degradation contribution value of the medium pore; and the fractal dimension of the large pore is multiplied by the attenuation correlation slope to obtain the degradation contribution value of the large pore. These three degradation contribution values ​​quantify the contribution of pore space complexity at different scale levels to the attenuation of bearing capacity.

[0126] It is understood that the overall load-bearing capacity degradation is the sum of the degradation contribution values ​​of small pores, medium pores, and large pores. This value comprehensively reflects the cumulative degradation effect of pores of all sizes on the material's load-bearing capacity under the combined action of fractal structure and stress concentration.

[0127] Preferably, the load-bearing performance degradation ratio is obtained by dividing the overall load-bearing performance degradation degree by the initial load-bearing capacity response value recorded when the material specimen is subjected to the same preset load before heating. The resulting quotient is the load-bearing performance degradation ratio. This ratio normalizes the degradation degree to a dimensionless proportional range based on the initial load-bearing capacity.

[0128] In one embodiment, the determination process of the multi-scale fractal influence coefficient is as follows: The fractal dimensions of small pores, medium pores, and large pores are extracted from the multi-scale fractal dimension; the three are added together and divided by three to obtain the mean fractal dimension; the load-bearing capacity degradation ratio is multiplied by the mean fractal dimension, and the resulting product is the multi-scale fractal influence coefficient of thermal damage on material toughness. Specifically, the calculation formula for the multi-scale fractal influence coefficient is as follows: K fractal =(ΔF total / F_0)×[(D s +D m +D l ) / 3], where K fractal ΔF is the multi-scale fractal influence coefficient. total The overall load-bearing capacity degradation is represented by F0, the initial load-bearing capacity response value is D. s D m D l These are the fractal dimensions of small, medium, and large pores, respectively. This influence coefficient simultaneously contains information on the proportion of load-bearing capacity degradation and the overall level of pore space complexity; a larger value indicates a stronger multi-scale fractal influence of thermally induced damage on the material's toughness.

[0129] Step S106: Using the multi-scale fractal influence coefficient as the damage weighting factor, a damage evolution path model is constructed in the digital twin platform. When the stress concentration point density in the high-risk area of ​​the fusion path continues to increase, the simulation boundary conditions are adjusted to obtain the final toughness assessment result.

[0130] Using the multi-scale fractal influence coefficient as a damage weighting factor, a three-dimensional finite element mesh of the material specimen is established in a digital twin platform. The three-dimensional spatial coordinates of all high-risk areas of the fusion path in the porosity fusion trend diagram are mapped to the corresponding elements of the finite element mesh. For each mapped mesh element, the preset elastic modulus of the element is divided by the damage weighting factor to obtain the reduced elastic modulus, resulting in a damage-weighted finite element mesh. Based on the damage-weighted finite element mesh, initial simulation boundary conditions are set. These conditions include temperature and mechanical loads applied to the outer surface of the mesh. The temperature load is taken from the temperature values ​​at each acquisition time in the surface temperature time series curve, and the mechanical load is taken from the preset load values ​​in the bearing capacity test. The initial simulation boundary conditions drive the finite element mesh to simultaneously solve for the temperature and stress fields, outputting the stress and displacement values ​​of each mesh element at each simulation time, thus obtaining the initial simulation results of damage evolution. Figure 3 As shown, based on the initial simulation results of the damage evolution, the simulated stress concentration point density of mesh elements in the high-risk area of ​​the fusion path is statistically analyzed at each simulation time. If the simulated stress concentration point density shows a continuous increasing trend for multiple consecutive simulation times and the increase exceeds a preset density growth threshold, the current simulation boundary conditions are adjusted. The adjustment method is to multiply the density increase by a preset load adjustment coefficient to obtain an additional temperature load increment, and then add the additional temperature load increment to the mesh element corresponding to that region. The adjusted boundary conditions are used to re-drive the finite element mesh solution to obtain updated damage evolution simulation results. Based on the updated damage evolution simulation results, the residual bearing capacity values ​​of all mesh elements at the final simulation time are extracted. The residual bearing capacity values ​​are divided by the initial bearing capacity response values ​​to obtain the residual bearing ratio. The residual bearing ratio is then divided by the damage weighting factor to obtain the final toughness assessment result.

[0131] Specifically, the method for determining and the range of values ​​for the load adjustment coefficient reflect the feedback intensity of stress concentration point density changes on the local temperature field. Its physical meaning is: when the stress concentration point density increases by 1 unit, the additional temperature load increment (unit: °C) required to be applied to that area. The load adjustment coefficient is related to the thermal conductivity and coefficient of thermal expansion of the material. Materials with lower thermal conductivity and higher coefficient of thermal expansion require a larger load adjustment coefficient to reflect the local heat accumulation effect. Specifically, the load adjustment coefficient for asphalt mixtures ranges from 2.0 to 5.0 °C·mm³ / piece, with a recommended value of 3.5 °C·mm³ / piece (based on asphalt's thermal conductivity of approximately 0.7 W / (m·K) and linear expansion coefficient of approximately 2 × 10⁻⁻⁻⁻⁶). 5 / ℃); Cement concrete: The load adjustment factor ranges from 0.5 to 2.0℃·mm³ / piece, with a recommended value of 1.2℃·mm³ / piece (based on the concrete thermal conductivity of approximately 1.5W / (m·K) and linear expansion coefficient of approximately 1×10⁻). 5 / ℃); Cement mortar: The load adjustment coefficient ranges from 0.8 to 2.5℃·mm³ / piece, with a recommended value of 1.5℃·mm³ / piece.

[0132] Specifically, the load adjustment coefficient can be calibrated through the following experiment: In the digital twin platform, the measured residual bearing capacity value is used as the target value, and the load adjustment coefficient is adjusted using an iterative optimization algorithm to minimize the error between the simulated residual bearing capacity value and the measured value. The load adjustment coefficient at this time is the calibration value of the material.

[0133] The density growth threshold is used to trigger dynamic adjustment of the simulated boundary conditions. This threshold is set to 10%-20% of the density of stress concentration points at the previous simulation time. For example, if the simulated stress concentration point density at the previous time is 0.05 points / mm³, then the density growth threshold is set to 0.005-0.01 points / mm³. When the density increment at three or more consecutive simulation times exceeds this threshold, boundary condition adjustment is triggered.

[0134] In one embodiment, the process of establishing the three-dimensional finite element mesh is as follows: In a digital twin platform, according to the actual geometric dimensions and shape of the material specimen, the three-dimensional space of the specimen is divided into several non-overlapping hexahedral mesh elements. Each mesh element has independent spatial coordinates and a preset elastic modulus value. The three-dimensional spatial coordinates of the high-risk areas of the fusion path in the porosity fusion trend diagram are matched with the element coordinates of the finite element mesh. The successfully matched mesh elements are the damage mapping elements.

[0135] Specifically, the elastic modulus reduction process is as follows: For each damaged mapping element, the original elastic modulus value of the element is divided by the damage weighting factor to obtain the reduced elastic modulus. When the damage weighting factor is greater than 1, the reduced elastic modulus decreases, correspondingly weakening the element's resistance to deformation in the simulation; when the damage weighting factor is less than 1, its reciprocal is taken before multiplying it by the original elastic modulus to ensure that the elastic modulus of the damaged region always decreases. The elastic modulus characterizes the material's ability to resist elastic deformation; dividing by a damage weighting factor greater than 1 reduces the value, correspondingly weakening the element's resistance to deformation in the simulation. Unmapped mesh elements retain their original elastic modulus, and the damage-weighted finite element mesh is formed after all elements are updated.

[0136] It should be noted that the temperature load and mechanical load in the initial simulated boundary conditions are both derived from the measured data obtained above: the temperature load is extracted from the temperature values ​​corresponding to each acquisition time from the surface temperature time series curve, and the mechanical load is taken from the preset load values ​​applied to the surface of the specimen during the bearing capacity test.

[0137] In one possible implementation, the process of simultaneously solving for the temperature and stress fields is as follows: At each simulation moment, based on the temperature load applied to the outer surface of the finite element mesh, the internal temperature values ​​of all mesh elements are solved one by one according to the heat conduction equation, obtaining the temperature field distribution of each element at that moment. Subsequently, based on the thermal expansion effect generated by the temperature field distribution, combined with the elastic modulus value of each element and the applied mechanical load, the stress and displacement values ​​of all elements are solved one by one according to the elastic equilibrium equation, obtaining the stress field distribution of each element at that moment.

[0138] For example, in the simulation of asphalt mixture specimens, when the external surface temperature load increases, the temperature of each element within the mesh is transferred layer by layer through heat conduction. The resulting thermal expansion causes mutual compression between the elements. After the superimposed mechanical load, the stress value borne by each element changes accordingly. The temperature and stress fields are solved sequentially for all simulation moments, outputting the stress and displacement values ​​of each element at each moment, forming the initial simulation results of damage evolution. Furthermore, the method for monitoring the simulated stress concentration point density is as follows: at each simulation moment, the number of mesh elements with stress values ​​exceeding a preset stress threshold in the high-risk area of ​​the fusion path is counted, and this number is divided by the volume of that area to obtain the simulated stress concentration point density at that moment.

[0139] For example, the triggering and execution process of the dynamic adjustment of the boundary conditions is as follows: The density of simulated stress concentration points at multiple consecutive simulation moments is arranged in chronological order, and the density values ​​of two adjacent simulation moments are compared moment by moment. If the density value at a later moment is continuously greater than that at a previous moment and the increase in density exceeds a preset density growth threshold, then the boundary condition adjustment is triggered. The specific adjustment method is as follows: the increase in density is multiplied by a preset load adjustment coefficient, and the resulting product is the additional temperature load increment. This additional temperature load increment is added to the outer surface of the mesh element corresponding to the region where the adjustment is triggered. The load adjustment coefficient is a fixed value pre-calibrated according to the material type and is used to control the magnitude of the temperature load increment. After the addition is completed, the temperature field and stress field are solved again on the finite element mesh with the adjusted boundary conditions.

[0140] It is understood that the result obtained by the re-solution is the updated damage evolution simulation result, which reflects the redistribution of stress and displacement inside the material after the boundary conditions are adjusted.

[0141] Preferably, the residual bearing capacity value is obtained by extracting the stress values ​​of all mesh elements at the final simulation time from the updated damage evolution simulation results, summing the stress values ​​of all elements and converting them into the mechanical load value that the whole specimen can withstand. This value is the residual bearing capacity value.

[0142] In one embodiment, the final toughness assessment result is determined by dividing the residual load-bearing capacity value by the initial load-bearing capacity response value to obtain the residual load-bearing ratio, and then dividing the residual load-bearing ratio by the damage weighting factor. The resulting value is the final toughness assessment result. The larger the value, the higher the level of toughness retained by the material after experiencing thermal damage.

[0143] Step S107: Based on the final toughness assessment results and the damage evolution path model, update the correlation database between heating rate, multi-scale fractal dimension and damage degree, and determine the optimized boundary conditions for subsequent thermal damage evolution simulation.

[0144] Based on the final toughness assessment result and the overall load-bearing capacity degradation degree, the current heating rate level, the multi-scale fractal dimension, the overall load-bearing capacity degradation degree, and the toughness assessment result are treated as a linked record and written into a pre-established association database. If a record with the same heating rate level already exists in the association database, the current linked record replaces the original record, completing the feedback update of the association database. According to the updated association database, the overall load-bearing capacity degradation degree and toughness assessment result corresponding to all heating rate levels are retrieved. The heating rate level with the largest overall load-bearing capacity degradation degree is matched, and the temperature load and load adjustment coefficient used in the simulation process for this level are extracted to determine the optimized boundary conditions for subsequent thermally induced damage evolution simulation.

[0145] In one implementation, the associated record consists of four fields: the heating rate level field records the constant heating rate value used in the current test; the multi-scale fractal dimension field records the combined values ​​of the fractal dimension of the small hole, the fractal dimension of the medium hole, and the fractal dimension of the large hole; the overall load-bearing capacity degradation degree field records the sum of the degradation contribution values ​​of the three scale levels; and the toughness assessment result field records the value of the residual load-bearing ratio divided by the damage weighting factor.

[0146] Specifically, the associated database uses the heating rate level as an index field, with each associated record corresponding to an independent heating rate level. During writing, the database is indexed according to the heating rate level in the current associated record. If a record with the same heating rate level already exists, the values ​​of the three fields—multi-scale fractal dimension, overall load-bearing capacity degradation degree, and toughness assessment result—of the original record are replaced with the corresponding values ​​in the current associated record. If no matching level is found, the current associated record is appended to the database as a new entry. Further, the process of determining the optimization boundary conditions is as follows: all associated records in the updated associated database are traversed, and the value of the overall load-bearing capacity degradation degree field in each record is read. The associated record with the largest value is located, and the heating rate level corresponding to this record represents the most severely damaged operating condition.

[0147] For example, from the simulation process corresponding to the most severe damage condition, the temperature load value and load adjustment coefficient value used at that level are extracted, and the two are used as the optimized boundary conditions for subsequent thermal damage evolution simulation, so that subsequent simulations prioritize the refinement and adjustment of boundary conditions for the heating rate condition with the highest degree of degradation.

[0148] If the technical solution of this application involves personal information, the product using this technical solution has clearly informed the user of the personal information processing rules and obtained the user's voluntary consent before processing the personal information. If the technical solution of this application involves sensitive personal information, the product using this technical solution has obtained the user's separate consent before processing the sensitive personal information, and also meets the requirement of "express consent". For example, at personal information collection devices such as cameras, clear and prominent signs are set up to inform users that they have entered the scope of personal information collection and that personal information will be collected. If an individual voluntarily enters the collection scope, it is deemed that they have agreed to the collection of their personal information; or on the personal information processing device, the personal information processing rules are clearly informed through signs / information, and authorization is obtained through pop-up information or by asking the individual to upload their personal information; wherein, the personal information processing rules may include information such as the personal information processor, the purpose of personal information processing, the processing method, and the types of personal information processed.

[0149] The above embodiments are merely one of the preferred embodiments of the present invention and should not be used to limit the scope of protection of the present invention. Any modifications or refinements made to the main design concept and spirit of the present invention that are not of substantial significance, but solve the same technical problem as the present invention, should be included within the scope of protection of the present invention.

Claims

1. A method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins, characterized in that, include: Real-time data on surface temperature changes of road materials at different heating rates are collected, and thermal stress concentration areas are identified by combining the data with preset temperature threshold ranges to obtain an initial thermal stress distribution map. Based on the initial thermal stress distribution map, a CT scan was used to perform a three-dimensional scan of the material interior to extract the spatial distribution density and position coordinates of small, medium and large pores, and to determine the multi-scale fractal dimension and pore size distribution density. The spacing between the two peaks of the pore distribution curve is analyzed by the multi-scale fractal dimension and the pore size distribution density. When the spacing between the two peaks is lower than a preset threshold, it is marked as a high-risk area of ​​the fusion path, and a pore fusion trend map is obtained. By using the pore fusion trend map, the abnormal decrease in the number of pores in the middle pore interval is extracted, and the situation where small pores directly fuse across the middle pores to form large pores is determined. Based on this, the damage evaluation benchmark is dynamically adjusted, and the distribution density of stress concentration points is determined. Based on the stress concentration point distribution density, loads were applied to the material specimens at each acquisition time, and the bearing capacity response values ​​were recorded. The difference in bearing capacity response values ​​between adjacent acquisition times is calculated and divided by the corresponding temperature difference to obtain the bearing capacity attenuation gradient, thus forming an attenuation gradient sequence. Extract the acquisition time of the abrupt change in the attenuation gradient value from the attenuation gradient sequence, and pair the stress concentration point distribution density at that time with the attenuation gradient value; Linear regression was performed with stress concentration point distribution density as independent variable and attenuation gradient value as dependent variable, and the slope of the fitted line was used as the attenuation correlation slope. The fractal dimensions of the small pores, medium pores, and large pores at the end of the heating process are extracted from the multi-scale fractal dimension and multiplied by the attenuation correlation slope to obtain the degradation contribution values ​​of the small pores, medium pores, and large pores. The overall load-bearing capacity degradation level is obtained by adding the three degradation contribution values; The overall load-bearing capacity degradation degree is divided by the initial load-bearing capacity response value before heating to obtain the load-bearing capacity degradation ratio. The load-bearing capacity degradation ratio is multiplied by the mean of the small hole fractal dimension, medium hole fractal dimension and large hole fractal dimension to determine the multi-scale fractal influence coefficient. Using the multi-scale fractal influence coefficient as a damage weighting factor, a damage evolution path model is constructed in the digital twin platform. When the stress concentration point density in the high-risk area of ​​the fusion path continues to increase, the simulation boundary conditions are adjusted to obtain the toughness assessment results. By using the toughness assessment results and the damage evolution path model, the correlation database between the heating rate, the multi-scale fractal dimension, and the degree of damage is updated to determine the optimal boundary conditions for subsequent thermal damage evolution simulation.

2. The method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins as described in claim 1, characterized in that, The obtained initial thermal stress distribution map includes: Infrared thermal imager was used to collect the surface temperature of road material specimens frame by frame. The surface temperature value of each frame was recorded according to the heating rate setting to generate a surface temperature time series curve. The temporal temperature difference between adjacent frames is calculated based on the surface temperature time-series curve to obtain a temporal temperature difference sequence. Extract continuous frame segments from the time-series temperature difference sequence where the time-series temperature difference value continuously exceeds a preset frame number threshold and is higher than the temperature threshold range, and mark the surface spatial location corresponding to the continuous frame segments as candidate regions for thermal stress concentration. For the candidate region of thermal stress concentration, temperature values ​​of sampling points are extracted along the horizontal and vertical directions, and the differences at each point are calculated to form gradient distribution characteristics; When the difference between horizontal and vertical points in the gradient distribution features both exceed the temperature threshold range, the corresponding candidate region is marked as a high thermal stress concentration region. Based on the spatial coordinates and gradient distribution characteristics of the high thermal stress concentration area, the average value of the point-by-point difference is taken on the two-dimensional coordinate plane, and the thermal stress intensity level is divided according to the average value and mapped to the corresponding position to obtain the initial thermal stress distribution map.

3. The method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins as described in claim 1, characterized in that, The determination of the multi-scale fractal dimension and aperture distribution density includes: Based on the spatial coordinates of the high thermal stress concentration area in the initial thermal stress distribution map, a layer-by-layer three-dimensional scan was performed using an industrial CT scanner to obtain voxel slices. Voxels with gray values ​​below the pore gray value threshold are marked, and their three-dimensional spatial coordinates and voxel volumes are extracted to form the original pore voxel set. For the original set of pore voxels, the equivalent diameter of the pores is calculated based on the number of voxels, and the pores are divided into small pore interval, medium pore interval and large pore interval according to the pore diameter threshold. Using voxel grid cells as statistical windows, the proportion of pores at each scale level is statistically analyzed window by window to obtain the pore window density ratio. The box counting method was used to process the pores at each scale level, the number of non-empty grids was counted step by step, and linear regression was performed. The absolute value of the slope of the fitted line was used as the fractal dimension of the corresponding scale, and the combination was used to obtain the multi-scale fractal dimension. The pore window density ratios at each scale level are aligned and summarized according to three-dimensional spatial coordinates to form the pore size distribution density.

4. The method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins as described in claim 1, characterized in that, The obtained pore fusion trend map includes: Based on the pore size distribution density, an overall pore distribution curve is plotted with the equivalent pore diameter as the horizontal axis and the pore frequency of the corresponding diameter interval as the vertical axis. Locate the local maxima of frequency on the pore distribution curve, mark the local maxima with smaller horizontal axis as small pore peaks, and mark the local maxima with larger horizontal axis as large pore peaks to obtain the bi-peak positioning coordinates; Calculate the difference in the equivalent pore diameter corresponding to the peak value of the small pore and the peak value of the large pore to obtain the original value of the double peak spacing; Extract the fractal dimension of the scale level to which the small hole peak belongs and the fractal dimension of the scale level to which the large hole peak belongs from the multi-scale fractal dimension, and divide the fractal dimension of the large hole peak to the fractal dimension of the small hole peak to obtain the fractal ratio value. Multiplying the fractal ratio by the original bimodal spacing value yields the fractal-corrected bimodal spacing. When the bimodal spacing after fractal correction is lower than the bimodal spacing threshold, the corresponding three-dimensional spatial location is marked as a high-risk area of ​​the fusion path; Threshold comparison and marking are performed on all spatial regions covered by the pore distribution curve one by one, and all high-risk areas of the fusion path are mapped according to three-dimensional spatial coordinates to obtain a pore fusion trend map.

5. The method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins as described in claim 1, characterized in that, Determining the stress concentration point distribution density includes: By using the spatial location of the high-risk area of ​​the fusion path in the pore fusion trend map, the current pore frequency of the middle pore interval in that area and the pore frequency of the middle pore interval at the previous acquisition time are extracted from the pore size distribution density, and the change in the middle pore frequency is calculated. Simultaneously, the current pore frequency of the large-pore interval in the high-risk area is extracted and compared with the pore frequency of the large-pore interval at the previous acquisition time, and the large-pore frequency increment is calculated. When the change in the frequency of the central hole is negative and the absolute value exceeds the threshold of the sudden drop, the corresponding area is marked as the abnormal decrease area of ​​the central hole. For the abnormal decrease region of the middle hole, the frequency increment of the large hole is compared with the absolute value of the frequency change of the middle hole. When the frequency increment of the large hole reaches a preset ratio threshold, it is determined that the small hole directly merges across the middle hole to form a large hole, and the region is marked as a cross-level fusion region. Within the cross-level fusion region, the damage evaluation benchmark is adjusted to the cross-level fusion extension of pores; Based on the spatial range of the cross-level fusion region, the three-dimensional spatial coordinates of each pore in the large-pore interval are located, and the coordinate position of each large-pore is marked as a stress concentration point. The total number of stress concentration points at each acquisition time within the high-risk area of ​​the fusion path is counted, and then divided by the volume of the spatial range to determine the stress concentration point distribution density sequence.

6. The method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins as described in claim 1, characterized in that, The obtained toughness assessment results include: The multi-scale fractal influence coefficient is used as the damage weighting factor to establish a three-dimensional finite element mesh in the digital twin platform; Map the three-dimensional spatial coordinates of the high-risk areas of the fusion path in the porosity fusion trend diagram to the corresponding elements of the finite element mesh; For the mapped mesh element, the preset elastic modulus is divided by the damage weighting factor to obtain the reduced elastic modulus, thus forming a damage-weighted finite element mesh. Initial simulation boundary conditions are set, including temperature load and mechanical load. The initial simulation boundary conditions are used to drive the finite element mesh to solve for the temperature field and stress field, and the initial simulation results of damage evolution are obtained. At each simulation moment, the density of simulated stress concentration points in the high-risk area of ​​the fusion path is statistically analyzed. When the density of simulated stress concentration points continues to increase for multiple consecutive simulation moments and the increase exceeds the density growth threshold, the simulation boundary conditions are adjusted, additional temperature load increments are added, and the finite element mesh is re-driven to obtain updated damage evolution simulation results. Based on the updated damage evolution simulation results, the residual bearing capacity values ​​of all grid elements at the final simulation time are extracted, the residual bearing capacity ratio is calculated, and the residual bearing capacity ratio is divided by the damage weight factor to obtain the toughness assessment results.

7. The method for simulating the thermally induced damage evolution and assessing the toughness of road materials by integrating multi-scale fractal features and digital twins as described in claim 1, characterized in that, The determination of the optimized boundary conditions for subsequent thermally induced damage evolution simulation includes: The current heating rate level, the multi-scale fractal dimension, the overall load-bearing capacity degradation degree, and the toughness assessment result are recorded as a single associated record and written into the associated database. When a record with the same heating rate already exists in the associated database, the existing record is replaced with the current associated record to complete the feedback update; Based on the updated associated database, retrieve the overall load-bearing capacity degradation and toughness assessment results corresponding to all heating rate levels; Match the heating rate setting with the highest overall load-bearing capacity degradation value, extract the temperature load and load adjustment coefficient used for the heating rate setting, and determine the optimized boundary conditions for subsequent thermal damage evolution simulation.

Images