Anti-rutting regulation and control method for high-mixing-amount milling material regenerated warm mix asphalt

Abstract

The invention discloses an anti-rut regulation and control method for regenerated warm mix asphalt of a large-mixing-amount milling material, and relates to the technical field of automatic control, and the method comprises the following steps: step 1, carrying out a needle penetration test and a softening point test on old asphalt in the milling material to obtain an old asphalt needle penetration value and an old asphalt softening point value; performing a melting point test on the to-be-doped warm mixing agent to obtain a melting point value of the warm mixing agent; 2, establishing a two-dimensional rectangular calculation domain, dividing the two-dimensional rectangular calculation domain into a first calculation section, a second calculation section and a third calculation section in the horizontal direction, and extracting anti-rutting performance evaluation indexes after iteration is completed; and 3, comparing the anti-rutting performance evaluation index with a preset target threshold value, adjusting the matching parameter according to a comparison result, returning to the step 2 for re-simulation until a target requirement is met, and determining a final matching parameter. According to the method, coupling simulation of diffusion evolution and high-temperature viscous flow behaviors of new and old asphalt interfaces is realized, and an optimal matching scheme can be quickly screened in a virtual environment.

Description

Technical Field

[0001] This invention relates to the field of automatic control technology, and in particular to a method for controlling rutting resistance in warm-mix asphalt with a large amount of milled material, which belongs to a type of industrial control software. Background Technology

[0002] In the development of recycled asphalt mixture technology, increasing the milled asphalt content has always been a key focus for researchers. Early recycling technologies, limited by production processes and mix design methods, typically controlled the milled asphalt content between 15% and 25%. In recent years, with continuous advancements in recycling technology, some projects have increased the milled asphalt content to 40% or even over 50%, maximizing the resource utilization of the milled asphalt. However, the old asphalt in the milled asphalt has undergone significant aging after long-term service, exhibiting reduced penetration, increased softening point, and decreased ductility, resulting in an overall hard and brittle characteristic. As the milled asphalt content increases, the proportion of aged asphalt in the mixture also increases, posing challenges to the high-temperature stability, low-temperature crack resistance, and fatigue durability of the recycled mixture.

[0003] Existing methods for designing recycled asphalt mixtures primarily rely on laboratory tests and empirical formulas. Designers first select an appropriate grade of new asphalt based on the aging degree of the old asphalt, then estimate the blending ratio of new and old asphalt using penetration or viscosity blending formulas. Finally, they prepare specimens for Marshall stability or rutting tests to verify the rationality of the design. This test-based design method has significant limitations. On the one hand, laboratory tests are time-consuming and costly, making it difficult to screen and optimize a large number of mix design schemes in a short period. On the other hand, traditional testing methods can only obtain macroscopic mechanical performance indicators, failing to reveal the microscopic fusion mechanism of the new and old asphalt interface and the rheological behavior of asphalt under load, thus lacking theoretical guidance in the design process. Summary of the Invention

[0004] The purpose of this invention is to provide a method for controlling the rutting resistance of recycled warm-mix asphalt with a large amount of milled material. It realizes the coupled simulation of the diffusion evolution of the new and old asphalt interface and the high-temperature viscous flow behavior, establishes a multi-dimensional evaluation system for rutting resistance performance, and can quickly screen and optimize the mix design in a virtual environment, which can significantly shorten the design cycle and reduce the test cost, providing scientific guidance for the mix design of recycled warm-mix asphalt mixtures with a large amount of milled material.

[0005] To solve the above-mentioned technical problems, the present invention provides the following technical solution: A method for controlling rutting resistance in warm-mix asphalt with high admixture of milled aggregate includes the following steps: Step 1: Conduct penetration and softening point tests on the old asphalt in the milled material to obtain the penetration value and softening point value of the old asphalt. Conduct penetration and softening point tests on the new asphalt to be added to obtain the penetration value and softening point value of the new asphalt. Conduct melting point tests on the warm mix additive to be added to obtain the melting point value of the warm mix additive. Step 2: Establish a two-dimensional rectangular computational domain and divide it horizontally into a first computational segment, a second computational segment, and a third computational segment. The first computational segment represents the old asphalt phase region, the third computational segment represents the new asphalt phase region, and the second computational segment represents the transition region between the old and new asphalt interfaces. Establish the phase field sequence parameter distribution to represent the phase state properties of asphalt, and establish a lattice Boltzmann velocity discrete lattice to simulate asphalt flow behavior. Perform alternating iterative calculations of the phase field update sub-step and the lattice Boltzmann collision migration sub-step, apply external load conditions simulating rutting formation, and extract anti-rutting performance evaluation indicators after iteration. Step 3: Compare the anti-rutting performance evaluation index with the preset target threshold, adjust the mix proportion parameters according to the comparison results, and return to Step 2 to re-simulate until the target requirements are met, and determine the final mix proportion parameters.

[0006] Furthermore, step one also includes: using a laser particle size analyzer to test the particle size distribution of aggregates in the milled material to obtain aggregate particle size distribution data; using an X-ray computed tomography (CT) scanner to scan the milled material sample to obtain a three-dimensional spatial distribution image of the internal pores of the milled material, and extracting the initial porosity value and the average pore diameter value of the milled material from the three-dimensional spatial distribution image.

[0007] Furthermore, in step two, the phase field sequence parameter ranges from zero to one in a continuous interval. When the phase field sequence parameter is zero, it indicates that the corresponding position is entirely old asphalt phase. When the phase field sequence parameter is one, it indicates that the corresponding position is entirely new asphalt phase. When the phase field sequence parameter is between zero and one, it indicates that the corresponding position is in a mixed transition state between old and new asphalt. When initializing the phase field sequence parameter, the initial value of the phase field sequence parameter at all positions in the first calculation section is set to zero, the initial value of the phase field sequence parameter at all positions in the third calculation section is set to one, and the initial value of the phase field sequence parameter at each position in the second calculation section is set to a transitional distribution that gradually changes from zero to one in the horizontal direction.

[0008] Furthermore, in step two, the lattice Boltzmann velocity discretization lattice adopts a D2Q9 lattice structure for velocity space discretization. The D2Q9 lattice structure sets nine discrete velocity directions at each lattice point. The nine discrete velocity directions include one zero velocity direction, four unit velocity directions along the positive and negative coordinate axes, and four unit velocity directions along the diagonal directions. Nine particle distribution functions are defined at each lattice point, and the nine particle distribution functions correspond to the nine discrete velocity directions. Each particle distribution function characterizes the number density of virtual particles moving along the corresponding discrete velocity direction.

[0009] Furthermore, step two also includes: setting a simulated temperature value that is higher than the softening point of the old asphalt, higher than the softening point of the new asphalt, and higher than the melting point of the warm mix agent; establishing a correlation rule between the phase field sequence parameter and the lattice Boltzmann local viscosity; determining the reference viscosity value of the old asphalt based on its penetration and softening point; determining the reference viscosity value of the new asphalt based on its penetration and softening point; and performing temperature correction on the reference viscosity values ​​of the old and new asphalt based on the simulated temperature value to obtain the temperature-corrected viscosity values ​​of the old and new asphalt. For any grid point within the computational domain, when the current phase field sequence parameter at the current grid point is zero, the local kinematic viscosity at the current grid point is set to the temperature-corrected old asphalt viscosity value; when the current phase field sequence parameter at the current grid point is one, the local kinematic viscosity at the current grid point is set to the temperature-corrected new asphalt viscosity value; when the current phase field sequence parameter at the current grid point is between zero and one, the local kinematic viscosity at the current grid point is set to the linear interpolation value between the temperature-corrected old asphalt viscosity value and the temperature-corrected new asphalt viscosity value, and the interpolation coefficient of the linear interpolation is the current phase field sequence parameter value.

[0010] Furthermore, the execution process of the phase field update sub-step in step two is as follows: For any target grid point within the computational domain, obtain the current phase field sequence parameter value of the target grid point, and obtain four adjacent phase field sequence parameter values: the current phase field sequence parameter values ​​of the target grid point in the positive horizontal direction, the current phase field sequence parameter values ​​of the adjacent grid points in the negative horizontal direction, the current phase field sequence parameter values ​​of the adjacent grid points in the positive vertical direction, and the current phase field sequence parameter values ​​of the adjacent grid point in the negative vertical direction. Calculate the arithmetic mean of the four adjacent phase field sequence parameter values ​​to obtain the adjacent average phase field sequence parameter value. Calculate the difference between the adjacent average phase field sequence parameter value and the current phase field sequence parameter value of the target grid point to obtain the phase field diffusion driving force. Obtain the current velocity vector at the target grid point, and calculate the product of the horizontal component of the current velocity vector and the difference between the current phase field sequence parameter values ​​of the adjacent grid points in the positive horizontal direction and the current phase field sequence parameter values ​​of the adjacent grid points in the negative horizontal direction. The horizontal convective transport is obtained. The vertical convective transport is obtained by multiplying the vertical component of the current velocity vector by the difference between the current phase field sequence parameter values ​​of adjacent grid points in the positive vertical direction and the current phase field sequence parameter values ​​of adjacent grid points in the negative vertical direction. The horizontal and vertical convective transport are added together to obtain the total convective transport. The diffusion driving force is multiplied by a preset phase field diffusion coefficient to obtain the diffusion update increment. The total convective transport is multiplied by a preset convective transport coefficient to obtain the convective update increment. The current phase field sequence parameter value of the target grid point is added to the diffusion update increment and then subtracted from the convective update increment to obtain the updated phase field sequence parameter value. When the updated phase field sequence parameter value is less than zero, it is corrected to zero. When the updated phase field sequence parameter value is greater than one, it is corrected to one.

[0011] Furthermore, the execution process of the lattice Boltzmann collision migration sub-step in step two is as follows: For any target lattice point within the computational domain, a collision operation is first performed. For each particle distribution function at the target lattice point, the local equilibrium distribution function value corresponding to the current particle distribution function is calculated based on the current macroscopic density value and the current velocity vector at the target lattice point. The difference between the current particle distribution function value and the local equilibrium distribution function value is calculated to obtain the deviation from equilibrium. The relaxation time value is calculated based on the local kinematic viscosity value at the target lattice point. The deviation from equilibrium is divided by the relaxation time value to obtain the collision adjustment amount. The current particle distribution function is then adjusted accordingly. Subtracting the collision adjustment amount from the value yields the particle distribution function value after the collision. After the collision operation is completed, a migration operation is performed, migrating the nine particle distribution functions after the collision at the target grid point to adjacent grid points along their respective discrete velocity directions. The particle distribution function after the collision corresponding to the zero velocity direction is retained at the target grid point. After the migration operation is completed, the particle distribution functions migrated to the current grid point are summarized for each grid point. The nine particle distribution function values ​​are added together to obtain the updated macroscopic density value. The nine particle distribution function values ​​are multiplied by their respective discrete velocity vectors, added together, and then divided by the updated macroscopic density value to obtain the updated flow velocity vector.

[0012] Furthermore, step two also includes setting boundary conditions at the boundaries of the two-dimensional rectangular computational domain: setting periodic boundary conditions at the upper and lower boundaries of the two-dimensional rectangular computational domain such that the particle distribution function migrating out from the upper boundary migrates into the lower boundary and the particle distribution function migrating out from the lower boundary migrates into the upper boundary; setting a constant flow velocity inlet boundary condition at the left boundary of the two-dimensional rectangular computational domain; and setting a constant pressure outlet boundary condition at the right boundary of the two-dimensional rectangular computational domain.

[0013] Furthermore, the method for applying the simulated rutting external load conditions in step two is as follows: the central section of the upper boundary of the two-dimensional rectangular computational domain is selected as the load application area, and a downward volumetric force source term is superimposed on the particle distribution function at each grid point within the load application area; the rutting resistance performance evaluation index includes the new and old asphalt interface fusion index, the simulated maximum rutting depth index, and the high-temperature flow activity index. The new and old asphalt interface fusion index is the percentage of grid points whose phase field sequence parameter values ​​are between the preset lower limit and upper limit of the mixing judgment value to the total number of grid points in the two-dimensional rectangular computational domain. The simulated maximum rutting depth index is the maximum value of the cumulative vertical displacement of the monitoring grid points directly below the load application area. The high-temperature flow activity index is the percentage of grid points whose velocity vector modulus is greater than the preset flow threshold to the total number of grid points in the two-dimensional rectangular computational domain.

[0014] Furthermore, the specific methods for adjusting the proportioning parameters in step three are as follows: when the fusion index of the new and old asphalt interface is less than the preset fusion target threshold, increase the amount of warm mix additive or increase the mixing temperature; when the simulated maximum rutting depth index is greater than the preset allowable rutting depth value, reduce the amount of milling material or increase the proportion of hard components in the new asphalt; when the high temperature flow activity index is greater than the preset upper limit of flow activity, increase the amount of polymer modifier in the new asphalt.

[0015] The method described in this invention has the following beneficial effects: First, this invention establishes a coupled simulation framework of the phase-field method and the lattice Boltzmann method, achieving a synergistic simulation of the diffusion evolution and high-temperature viscous flow behavior at the interface between new and old asphalt. The phase-field method uses continuously varying order parameters to describe the phase distribution of new and old asphalt, transforming the traditional sharp interface into a diffusion transition region with a certain thickness. This avoids the numerical difficulties of interface tracking and naturally simulates the mutual penetration and fusion process of new and old asphalt molecules under warm-mix conditions. The lattice Boltzmann method simulates the macroscopic flow behavior of asphalt through the collision migration rules of the particle distribution function. Its local computational characteristics make the algorithm easy to parallelize, resulting in significantly higher computational efficiency than traditional finite element or finite volume methods. The coupling of the two methods is achieved through the correlation between the phase-field order parameters and local kinematic viscosity, allowing the viscosity of the interface region to smoothly transition with phase changes, realistically reflecting the rheological characteristics of the new and old asphalt mixing region.

[0016] Secondly, this invention achieves a quantitative correlation between mix proportioning parameters and rutting resistance. By comparing simulation results with target thresholds and adjusting the content of milled aggregate, warm mix additive, new asphalt, and the proportion of new asphalt components accordingly, an iterative optimization closed loop based on simulation prediction is formed. This method overcomes the limitations of traditional mix design relying on numerous indoor tests, enabling rapid screening and optimization of mix proportioning schemes in a computer virtual environment, significantly shortening the design cycle and reducing testing costs. Simultaneously, the simulation process can intuitively demonstrate the evolution of the interface between new and old asphalt and the development law of rutting deformation, providing theoretical support for understanding the rutting resistance mechanism of recycled warm mix asphalt mixtures, transforming mix design from empirical trial and error to theoretical guidance. Attached Figure Description

[0017] Figure 1 Viscosity-temperature characteristic curves of different asphalt components and recycled mixtures, and schematic diagram of the action mechanism of warm mix agent provided in the embodiments of the present invention; Figure 2 A schematic diagram of the distribution of phase field sequence parameters and velocity vector field at a certain moment, obtained by coupled simulation based on the phase field method and the lattice Boltzmann method, as provided in an embodiment of the present invention; Figure 3 The biaxial curves showing the evolution of the anti-rutting performance evaluation index and interface fusion degree with the number of simulation iterations provided in the embodiments of the present invention. Detailed Implementation

[0018] A method for controlling rutting resistance in warm-mix asphalt with high admixture of milled aggregate includes the following steps: Step 1: Conduct penetration and softening point tests on the old asphalt in the milled material to obtain the penetration value and softening point value of the old asphalt. Conduct penetration and softening point tests on the new asphalt to be added to obtain the penetration value and softening point value of the new asphalt. Conduct melting point tests on the warm mix additive to be added to obtain the melting point value of the warm mix additive. Step 2: Establish a two-dimensional rectangular computational domain and divide it horizontally into a first computational segment, a second computational segment, and a third computational segment. The first computational segment represents the old asphalt phase region, the third computational segment represents the new asphalt phase region, and the second computational segment represents the transition region between the old and new asphalt interfaces. Establish the phase field sequence parameter distribution to represent the phase state properties of asphalt, and establish a lattice Boltzmann velocity discrete lattice to simulate asphalt flow behavior. Perform alternating iterative calculations of the phase field update sub-step and the lattice Boltzmann collision migration sub-step, apply external load conditions simulating rutting formation, and extract anti-rutting performance evaluation indicators after iteration. Step 3: Compare the anti-rutting performance evaluation index with the preset target threshold, adjust the mix proportion parameters according to the comparison results, and return to Step 2 to re-simulate until the target requirements are met, and determine the final mix proportion parameters.

[0019] Before obtaining physical property parameters, it is necessary to separate the old asphalt from the milled material recovered on-site for subsequent testing. Milled material typically consists of old asphalt, aggregate particles, and a small amount of impurities, all tightly bound together. Direct testing cannot accurately reflect the true performance state of the old asphalt. The trichloroethylene solvent extraction method can effectively separate the old asphalt from the aggregate surface. In practice, 200-300 grams of the milled material sample are weighed and placed in the filter paper tube of the extractor. 500 ml of trichloroethylene solvent is added to the extraction flask. The heating device is activated to boil the solvent and generate rising steam. The steam condenses in the condenser and drips into the filter paper tube, soaking the milled material sample. After the solvent dissolves the old asphalt, it flows back into the extraction flask through the bottom of the filter paper tube. This extraction process is repeated until the reflux liquid is colorless and transparent, indicating that the old asphalt has been fully dissolved and extracted. The asphalt solution in the extraction flask is transferred to a distillation apparatus and distilled at 160 degrees Celsius to remove the trichloroethylene solvent. The residue is the recovered old asphalt. Because trichloroethylene is somewhat toxic, operating it in a fume hood can effectively protect the safety of laboratory personnel. As an alternative, n-butanol can be used as the extraction solvent. While n-butanol has a slightly lower solubility for asphalt than trichloroethylene, it is less toxic; the extraction time can be extended to 1.5 times the original time to achieve the same effect.

[0020] The penetration test characterizes the hardness of asphalt under specified temperature conditions, and its value directly reflects the asphalt's resistance to external force intrusion. Recycled asphalt is heated to a fluid state and then poured into a standard penetration test dish. The dish has an inner diameter of 55 mm and a depth of 35 mm. During pouring, the asphalt level should be approximately 5 mm above the edge of the dish to compensate for cooling shrinkage. After pouring, the asphalt is cooled at room temperature for 1.5 hours to allow it to fully solidify. Then, the test dish, along with the asphalt sample, is placed in a constant-temperature water bath. The water bath temperature is strictly controlled at 25 degrees Celsius with a fluctuation range not exceeding ±0.1 degrees Celsius, and the constant-temperature time is no less than 1.5 hours to ensure uniform internal temperature of the sample. The test dish is then removed and placed on the test stage of the penetration tester. The standard needle is adjusted so that its tip just touches the surface of the asphalt sample. The standard needle has a mass of 100 grams and a cone angle of 9 degrees 12 minutes. The needle is released so that it penetrates vertically into the asphalt sample under its own weight for 5 seconds. The penetration depth is recorded as the old asphalt penetration value, in units of 0.1 mm. Three different locations on the same sample were selected for repeated testing. The distance between the three test points was no less than 10 mm, and the distance from the edge of the sample dish was no less than 10 mm. The arithmetic mean of the three test results was taken as the final penetration value of the old asphalt. Due to the volatilization of lightweight components and the aggregation of asphaltenes, the penetration value of old asphalt that has undergone long-term aging is usually between 15 and 40, while the penetration value of new asphalt that has not been aged is generally between 60 and 80, showing a significant difference. This difference is the experimental basis for setting different viscosity parameters for new and old asphalt in the subsequent phase-field method simulation.

[0021] refer to Figure 1 The horizontal axis in the graph represents the temperature variable. The unit is Celsius, using a linear scale; the vertical axis represents the dynamic viscosity variable. The unit is Pascal-second, using a double logarithmic scale to visually represent the exponential decay of asphalt viscosity with temperature. The figure includes three main characteristic curves, corresponding to old asphalt, new asphalt, and recycled asphalt with warm mix additives. The upper dashed line represents the viscosity-temperature curve of old asphalt. Due to long-term photo-oxidative and thermo-oxidative aging, the asphaltenes content in old asphalt increases significantly, and the light components volatilize, resulting in a significantly higher viscosity at the same temperature compared to other components, exhibiting higher stiffness and poorer flowability. The middle dotted line represents the viscosity-temperature curve of new asphalt added as a recycling agent. It is rich in aromatic and saturated components, has a relatively low viscosity, and mainly plays a softening and blending role in the system. The lower solid line represents the viscosity-temperature curve of the recycled asphalt mixture after adding the warm mix additive. This curve exhibits nonlinear rheological behavior within a specific temperature range. Specifically, a vertical dashed line in the figure corresponds to the melting point temperature of the warm mix additive. In this embodiment Set to 94 degrees Celsius.

[0022] At temperature Below At that time, the warm-mixing agent exists in solid form and has a limited effect on reducing the viscosity of the system; when the temperature... Rise and exceed At this time, the warm mix agent undergoes a phase change from solid to liquid, expanding in volume and producing a lubricating effect. This causes a significant step drop in the viscosity curve of the recycled asphalt, as shown in the figure. The inflection point and the viscosity decrease indicated by the arrow demonstrate this. This viscosity abrupt change mechanism is the core of warm mix technology, enabling the mixture to achieve sufficient workability at lower temperatures. The figure also shows another vertical dashed line corresponding to the simulated mixing temperature. In this embodiment, the temperature is set to 135 degrees Celsius. At this temperature, Significantly higher than This ensures that the warm mix additive completely melts to exert its viscosity-reducing effect, while also maintaining a temperature higher than the softening points of both old and new asphalt. This ensures that all components are in a Newtonian or near-Newtonian fluid state, satisfying the fundamental assumptions of fluid flowability in the coupled simulation of the phase-field method and the lattice Boltzmann method. Through comparison... The ordinate values ​​of the three curves can clearly quantify the effect of the warm mix additive on the viscosity reduction of aged asphalt, providing an accurate physical basis for setting the local kinematic viscosity parameters of each phase in subsequent numerical simulations.

[0023] The softening point test was conducted using the ring and ball method, which characterizes the high-temperature stability of asphalt by measuring the temperature at which it reaches a specified degree of deformation during the softening process. Old asphalt was heated to a fluid state and poured into a standard brass sample ring with an inner diameter of 15.9 mm and a height of 6.4 mm. After pouring, the surface was smoothed and cooled to room temperature. The sample ring was fixed to the support of the ring and ball apparatus, and a standard steel ball with a diameter of 9.53 mm and a mass of 3.5 g was placed in the center of the sample ring. The entire apparatus was immersed in a beaker containing distilled water, with a ring-shaped metal support at the bottom of the beaker to ensure that the sample ring was 50 mm below the water surface. The heating device was started, and the temperature was uniformly increased at a rate of 5 degrees Celsius per minute. As the temperature increased, the asphalt gradually softened, and the steel ball, under its own gravity, dragged the asphalt downwards. When the asphalt sank to contact the metal base plate below, the water bath temperature at this moment was recorded as the softening point value of the old asphalt. Two parallel tests were conducted, and the difference between the two test results should not exceed 1 degree Celsius. The average value was taken as the final softening point value of the old asphalt. The softening point of heavily aged asphalt typically ranges from 55 to 70 degrees Celsius, which is 5 to 15 degrees Celsius higher than that of new asphalt. A higher softening point means that the asphalt has a stronger ability to maintain its shape under high-temperature conditions, but it also indicates that the asphalt becomes more brittle and its low-temperature crack resistance decreases. In subsequent coupled simulations using the phase-field method and the lattice Boltzmann method, the softening point value is used to determine the lower limit of the simulation temperature. Only when the simulation temperature is higher than the softening point values ​​of both new and old asphalt can the asphalt exhibit a flow state, thus allowing interfacial diffusion and fusion to occur.

[0024] The new asphalt to be added is subjected to penetration and softening point tests using the same methods as the old asphalt. As a major component of the recycling agent, the performance indicators of the new asphalt directly affect the final road performance of the recycled mixture. In practical engineering, the selection of new asphalt needs to comprehensively consider the aging degree of the old asphalt. When the old asphalt penetration value is below 25 and the softening point value is above 60 degrees Celsius, it indicates severe aging, and soft new asphalt with a penetration value between 80 and 100 should be selected to compensate for the hardening of the old asphalt. When the old asphalt penetration value is between 25 and 40 and the softening point value is between 50 and 60 degrees Celsius, it indicates moderate aging, and ordinary new asphalt with a penetration value between 60 and 80 is sufficient. This method of selecting matching new asphalt based on the aging degree of the old asphalt ensures that the performance of the recycled asphalt mixture is close to the target design value. In this embodiment, the obtained penetration and softening point values ​​of the new asphalt, together with those of the old asphalt, are used to calculate the reference viscosity values ​​of the new and old asphalt in the phase-field method simulation. There is an empirical correlation between asphalt viscosity and penetration; a higher penetration value indicates softer asphalt and a lower corresponding viscosity value.

[0025] The melting point of the warm mix agent was tested using differential scanning calorimetry (DSC), a method that accurately captures the heat changes during the phase transition of the warm mix agent. A sample of 5 to 10 mg of warm mix agent was weighed and placed in an aluminum sample crucible. An empty aluminum crucible was used as a reference sample. Both crucibles were simultaneously placed in the sample chamber of the DSC. The temperature was increased from room temperature to 150°C at a rate of 10°C per minute. The instrument automatically recorded the heat flow difference between the sample and the reference sample as a function of temperature. When the warm mix agent undergoes a solid-to-liquid phase transition, it needs to absorb additional heat to overcome intermolecular forces. This process is represented by a downward-sloping endothermic peak on the heat flow curve. The temperature corresponding to the onset of the endothermic peak is the melting point of the warm mix agent. The melting point of commercially available warm mix agents is typically between 80°C and 120°C. In this embodiment, a warm mix agent with a melting point in the range of 90°C to 100°C is preferred. If the melting point of the warm-mix agent is too high, it will be difficult for it to fully melt and disperse during mixing, reducing the warm-mixing effect; if the melting point is too low, it will melt and clump prematurely during storage and transportation, affecting the metering accuracy. In the subsequent coupled simulation of the phase-field method and the lattice Boltzmann method, the simulated temperature must be higher than the melting point of the warm-mix agent to ensure that the warm-mix agent has completely melted under the simulated mixing conditions and exerts its viscosity-reducing effect.

[0026] As an optional implementation, the penetration test can also be conducted at 15 degrees Celsius or 30 degrees Celsius to obtain penetration values ​​at different temperatures. By measuring the penetration values ​​at multiple temperature points, the penetration index can be calculated, thus providing a more comprehensive evaluation of the temperature sensitivity of asphalt. Penetration Index The calculation involves a linear regression of the logarithm of the penetration value at different temperatures with temperature. A higher value indicates lower temperature sensitivity of the asphalt, and smaller performance fluctuations with temperature changes. For recycled warm-mix asphalt mixtures with high rutting resistance requirements, this value is preferred. Asphalt with a softening point value greater than -0.5. In addition, a glycerol bath can be used instead of a water bath for softening point testing. Glycerol has a boiling point of approximately 290 degrees Celsius and is suitable for testing modified asphalt with a softening point value exceeding 80 degrees Celsius. When the old asphalt is styrene-butadiene-styrene block copolymer modified asphalt, its softening point value may reach 75 to 90 degrees Celsius, in which case a glycerol bath must be used for testing. Optional methods for testing the melting point of warm mix additives include the direct observation method using a melting point apparatus. The warm mix additive powder is filled into a capillary tube and inserted into a heating furnace. The temperature is slowly increased at a rate of 1 degree Celsius per minute, and the melting point value is observed through a magnifying glass at the temperature at which the powder begins to melt. This method has simple equipment but its accuracy is slightly lower than that of differential scanning calorimetry (DSC).

[0027] The establishment of a two-dimensional rectangular computational domain is the foundation of the entire coupled simulation. A two-dimensional array space is allocated in computer memory to store the state information of each grid point within the computational domain. The horizontal length of the computational domain is set to 500 grid points, and the vertical height is set to 200 grid points. Each grid point represents a physical dimension of 1 micrometer. Therefore, the physical region corresponding to the entire computational domain is a rectangular area of ​​500 micrometers horizontally and 200 micrometers vertically. This scale is chosen based on the typical characteristics of the microstructure of asphalt mixtures. The width of the interface transition zone between new and old asphalt is typically between 50 and 150 micrometers, and the horizontal length of 500 micrometers can completely encompass the interface transition zone and the pure phase regions on both sides. When the computational domain is divided into three adjacent segments along the horizontal direction, the first computational segment occupies the range from the 1st to the 150th grid point, the second computational segment occupies the range from the 151st to the 350th grid point, and the third computational segment occupies the range from the 351st to the 500th grid point. The first calculation section characterizes the spatial area occupied by the old asphalt phase. In recycled warm-mix asphalt mixtures, the old asphalt originates from the aged asphalt film covering the surface of the milled material. Its molecular structure undergoes significant changes due to long-term photo-oxidative and thermo-oxidative aging, resulting in increased asphaltene content and decreased aromatic and saturated content, exhibiting overall high viscosity and high brittleness. The third calculation section characterizes the spatial area occupied by the new asphalt phase. As the main component of the recycling agent, the new asphalt is rich in lightweight components, capable of penetrating into the interior of the old asphalt and dissolving some aging products, restoring the fluidity and viscoelasticity of the old asphalt. The second calculation section, located between the old and new asphalt phases, characterizes the interfacial transition region formed after the two phases come into contact. In this region, the molecules of the old and new asphalt diffuse and permeate each other, which is a key area determining the recycling effect.

[0028] The introduction of the phase field order parameter allows the interface between new and old asphalt to be described in a continuously varying manner, avoiding the numerical difficulties encountered by traditional sharp interface methods when dealing with interface movement and deformation. The phase field order parameter is defined as follows: As a scalar field characterizing the phase properties of asphalt at any location within the computational domain, The value range of is a continuous interval from 0 to 1. When When the value is 0, it indicates that the corresponding position is completely occupied by the old bitumen phase; when... When the value is 1, it indicates that the corresponding position is completely occupied by the new bitumen phase; when A value between 0 and 1 indicates that the corresponding location is in a mixed transition state between old and new asphalt. The specific numerical value reflects the volume fraction of the new bitumen phase at that location. The advantage of this continuous processing method is that the interface is no longer a mathematical surface without thickness, but a transition region with a certain width, allowing physical quantities within the interface to change smoothly, thus ensuring the stability of the numerical calculation. When initializing the phase field order parameters, all grid points within the first calculation segment are traversed, and the values ​​at each grid point are... The initial value is set to 0; all grid points in the third calculation segment are traversed, and the initial value at each grid point is set to 1; for grid points in the second calculation segment, a linear interpolation method is used to set the initial value. The initial value, specifically, is located in the horizontal direction. At each grid point The initial value is calculated by... Subtracting 151 and dividing by 199 results in The value smoothly transitions from 0 at the 151st grid point to 1 at the 350th grid point.

[0029] The lattice Boltzmann method uses mesoscopic-scale particle distribution functions to describe the motion of fluids, simulating complex flow behaviors through simple collision and migration rules. This embodiment uses a D2Q9 lattice structure for velocity space discretization, where D2Q9 indicates that each lattice point in two-dimensional space has nine discrete velocity directions. The nine discrete velocity directions are specifically defined as follows: the 0th direction is the zero-velocity direction, corresponding to stationary particles; the 1st to 4th directions are the positive horizontal direction, positive vertical direction, negative horizontal direction, and negative vertical direction, respectively, corresponding to particles moving along the coordinate axes, with a velocity magnitude of 1 lattice unit per time step; the 5th to 8th directions are the upper right diagonal direction, upper left diagonal direction, lower left diagonal direction, and lower right diagonal direction, respectively, corresponding to particles moving along the diagonal, with a velocity magnitude of √2 lattice units per time step. Nine particle distribution functions are defined at each lattice point. to Each particle distribution function characterizes the number density of virtual particles moving along the corresponding discrete velocity direction. The advantage of using particle distribution functions instead of directly solving the macroscopic fluid equations is that the evolution rules of the particle distribution functions are local, and the update of each grid point depends only on its own and the information of its neighboring grid points. This characteristic makes the lattice Boltzmann method very suitable for parallel computing, and can make full use of the massive parallel computing capabilities of modern graphics processors.

[0030] Establishing a correlation rule between phase-field sequence parameters and lattice Boltzmann local viscosity is key to coupling the two methods. Based on the penetration and softening point values ​​of the old asphalt obtained in step one, the reference viscosity value of the old asphalt is calculated using empirical formulas. The reference viscosity of old asphalt with a penetration value of 25 and a softening point of 60 degrees Celsius is approximately 800 Pascals per second. Similarly, the reference viscosity of new asphalt is calculated based on its penetration and softening point values. A new asphalt with a penetration value of 70 and a softening point of 48 degrees Celsius corresponds to a reference viscosity of approximately 200 Pascals per second. (Set the simulation temperature value.) In this embodiment, the simulated temperature is set to 135 degrees Celsius. This temperature is higher than the softening point of the old asphalt, higher than the softening point of the new asphalt, and higher than the melting point of the warm mix additive, ensuring that both the old and new asphalt are in a flowable state under simulated conditions. The reference viscosity value is corrected based on the simulated temperature. The asphalt viscosity decreases significantly with increasing temperature. After temperature correction, the old asphalt viscosity is approximately 50 Pascals per second, and the new asphalt viscosity is approximately 8 Pascals per second. For any grid point within the computational domain, the phase field sequence parameter value at the current grid point is read. ,when When the value is 0, the local kinematic viscosity at the current grid point is... Set to the temperature-corrected viscosity value of the old asphalt, when When the value is 1, the local kinematic viscosity is set to the temperature-corrected new asphalt viscosity value. When the viscosity is between 0 and 1, linear interpolation is used to calculate the local kinematic viscosity, and the interpolation coefficient is... This viscosity interpolation method based on phase field order parameters can naturally reflect the gradual viscosity variation characteristics within the interfacial region, avoiding numerical oscillations caused by abrupt viscosity changes.

[0031] The phase field update sub-step is used to simulate the diffusion evolution process at the interface between new and old asphalt. For any target grid point in the computational domain, the current phase field sequence parameter value of the target grid point is first obtained. Then, the current phase field sequence parameter values ​​of the target grid point's adjacent grid points in the positive horizontal direction, negative horizontal direction, positive vertical direction, and negative vertical direction are obtained, resulting in four adjacent phase field sequence parameter values. The arithmetic mean of the four adjacent phase field sequence parameter values ​​is calculated to obtain the adjacent average phase field sequence parameter value. The difference between the adjacent average phase field sequence parameter value and the current phase field sequence parameter value of the target grid point is calculated to obtain the phase field diffusion driving force. This difference reflects the degree of inhomogeneity of phase properties between the target grid point and its surrounding environment. A positive difference indicates that the concentration of new asphalt in the surrounding area is higher than that at the target grid point, and the diffusion trend causes the new asphalt to penetrate towards the target grid point. The current velocity vector at the target grid point, calculated by the lattice Boltzmann collision migration sub-step, is obtained. The current velocity vector includes a horizontal component. and vertical component When calculating horizontal convective transport volume, The phase sequence parameter is calculated as the product of the current phase sequence parameter value at the grid point adjacent in the positive horizontal direction and the difference between the current phase sequence parameter value at the grid point adjacent in the negative horizontal direction. This calculation uses a central difference scheme to ensure second-order accuracy. Similarly, the vertical convective transport is calculated, and the two are added together to obtain the total convective transport. Convective transport describes the transport of phase sequence parameters caused by the macroscopic motion of the fluid; when the fluid travels from a high altitude... Regional flow to low When the phase field diffusion reaches a certain area, it will carry new asphalt components with it. The diffusion renewal increment is obtained by multiplying the phase field diffusion driving amount by a preset phase field diffusion coefficient. The time step was set to 0.001 square grid units, a value calculated based on the interdiffusion coefficient of asphalt components at high temperatures. The total convective transport was multiplied by a preset convective transport coefficient (set to 0.5) to obtain the convective update increment. The updated phase sequence parameter value was obtained by adding the diffusion update increment to the current phase sequence parameter value at the target grid point and then subtracting the convective update increment. This was done to prevent numerical errors from causing... If the value exceeds the physical reasonable range, the updated phase field sequence parameter value will be corrected to 0 if it is less than 0, and corrected to 1 if it is greater than 1.

[0032] The lattice Boltzmann collision migration sub-step is used to simulate the flow behavior of asphalt under external forces. The collision operation describes the interaction between particles at the lattice point, causing the particle distribution function to relax towards a local equilibrium state. For any target lattice point in the computational domain, collision operations are performed for each of the nine particle distribution functions. First, the current macroscopic density value is calculated based on the values ​​of the nine particle distribution functions at the target lattice point. The current flow velocity vector and the macroscopic density value are calculated by summing the values ​​of the nine particle distribution functions. Each component of the flow velocity vector is equal to the sum of the nine particle distribution function values ​​multiplied by their corresponding discrete velocity components, then divided by the macroscopic density value. Based on the current macroscopic density value and the current flow velocity vector, the local equilibrium distribution function value corresponding to each particle distribution function is calculated. The local equilibrium distribution function reflects the thermal equilibrium state of the particle distribution function under given macroscopic density and flow velocity conditions. The deviation from equilibrium is calculated by measuring the difference between the current particle distribution function value and the local equilibrium distribution function value; a larger deviation indicates a greater difference between the current state and the equilibrium state. The local kinematic viscosity value is obtained by querying the current phase field sequence parameter value at the target grid point, and the relaxation time value is calculated based on the local kinematic viscosity value. The relaxation time value is directly proportional to the local kinematic viscosity value. Higher viscosity means stronger internal fluid friction, a slower relaxation rate of the particle distribution function towards equilibrium, and a larger corresponding relaxation time value. Dividing the deviation from equilibrium by the relaxation time value yields the collision adjustment amount, and subtracting the collision adjustment amount from the current particle distribution function value gives the post-collision particle distribution function value. This single-relaxation-time collision model is simple and efficient, significantly reducing computational complexity while maintaining accuracy.

[0033] The migration operation propagates the particle distribution functions after a collision to adjacent grid points along their respective discrete velocity directions. Specifically, it migrates the particle distribution functions at the target grid point corresponding to the discrete velocities along the positive horizontal direction to adjacent grid points in the positive horizontal direction, the negative horizontal direction, the positive vertical direction, and the negative vertical direction. It also migrates the particle distribution functions corresponding to the discrete velocities along the four diagonal directions to their respective four diagonal adjacent grid points. The particle distribution function corresponding to the zero-velocity direction is retained at the target grid point. After the migration operation, for each grid point, the particle distribution functions migrated from all adjacent grid points to the current grid point, along with the zero-velocity direction particle distribution function retained at the current grid point, are summed. The updated macroscopic density value is obtained by adding the nine particle distribution function values ​​together. The updated flow velocity vector is obtained by multiplying each of the nine particle distribution function values ​​by its corresponding discrete velocity vector, summing them, and then dividing by the updated macroscopic density value.

[0034] refer to Figure 2 This figure presents the local physical state of a two-dimensional rectangular computational domain, revealing the microscopic mechanisms of the evolution of the interface between old and new asphalt and rutting formation. The variations in background color represent phase-field order parameters. The spatial distribution. According to the color mapping scale, the dark area (or blue area) corresponds to... A value close to 0 indicates a region of old asphalt phase aggregation, with high local viscosity; light-colored areas (or red areas) correspond to... A value close to 1 indicates a region of new asphalt phase aggregation with low local viscosity. Between these two pure phase regions lies a color-gradient transition zone, the interface diffusion layer between the old and new asphalt. Values ​​between 0 and 1 indicate that molecular-scale interpenetration and fusion have occurred between the old and new asphalt at this location. The figure not only displays the scalar field... The distribution is also superimposed with an array of white arrows to represent the velocity vector of the fluid particles. .

[0035] The direction of the arrow indicates the instantaneous direction of fluid motion, and the length of the arrow is related to the magnitude of the flow velocity vector. Proportional. As can be observed from the figure, directly below the region above the computational domain where an external load is applied, the velocity vector exhibits a significant downward component. This indicates that the asphalt material underwent vertical compression deformation under wheel load pressure; while on both sides of the load area, the velocity vector gradually turned to the horizontal direction, exhibiting an outward horizontal component. This indicates that the asphalt material undergoes shear flow and heaving on both sides after being compressed, which is a typical rheological characteristic of rutting distress in road surfaces. In particular, in the interfacial transition zone, the distribution of the velocity vector is significantly affected by the local viscosity gradient. Due to the lower viscosity of the nascent asphalt side, the velocity vector modulus on this side is generally greater than that on the unglued asphalt side, leading to non-uniform convective deformation at the interface. The bending and twisting morphology of the interface profile in the figure is precisely due to the convection term in the phase field evolution equation. Driven by this, real-time monitoring of the flow field and phase field allows for the calculation of the high-temperature flow activity index, i.e., the proportion of grid points with a statistical velocity vector modulus greater than a preset flow threshold. This quantitatively assesses the ability of asphalt mixtures with the current mix proportions to resist high-temperature shear deformation. This visualized simulation result intuitively reveals the intrinsic relationship between component distribution, viscosity differences, and macroscopic deformation.

[0036] Appropriate boundary conditions are set at the boundaries of the two-dimensional rectangular computational domain to simulate actual working conditions. Periodic boundary conditions are set at the upper and lower boundaries. The particle distribution function migrating out from the upper boundary migrates into the lower boundary, and vice versa. The periodic boundary conditions eliminate the obstruction effect of the upper and lower boundaries on the flow, allowing the simulation results to reflect the flow characteristics in an infinitely large area. A constant velocity inlet boundary condition is set at the left boundary. The fluid entering the computational domain from the left boundary has a preset inlet velocity value, set to 0.01 grid units per time step, which corresponds to the relative flow velocity of asphalt during actual mixing. A constant pressure outlet boundary condition is set at the right boundary, allowing the fluid to flow freely out of the computational domain without reflection.

[0037] When applying external load conditions simulating rutting formation, the central section from the 200th to the 300th grid point in the horizontal direction at the upper boundary is selected as the load application area. This section is 100 grid points wide, corresponding to a physical size of 100 micrometers, simulating the local action of the tire-road contact area. At each grid point within the load application area, a downward volumetric force term is superimposed onto the particle distribution function. The magnitude of the volumetric force term is determined based on a preset simulated wheel load pressure value. In this embodiment, the simulated wheel load pressure value is set to 0.7 MPa, which translates to approximately 0.0001 grid units for the volumetric force term. The application of volumetric force causes the asphalt below the load application area to experience continuous downward pressure. Under high temperature conditions, the asphalt undergoes viscous flow, gradually squeezing out to both sides and forming a depression below the load area. This process is consistent with the actual rutting formation mechanism of road surfaces.

[0038] The phase field update sub-step and the lattice Boltzmann collision migration sub-step are repeated until the preset total number of iterations is reached. Determining the total number of iterations requires establishing a conversion relationship between lattice time steps and actual physical time. Based on the similarity principle, one lattice time step corresponds to approximately 0.001 seconds of actual physical time. The loading time in the actual rutting test is 60 minutes, or 3600 seconds, resulting in a total of 3,600,000 iterations. Considering computational efficiency, a time step amplification technique can be used to reduce the total number of iterations to 36,000, with each lattice time step corresponding to 0.1 seconds of actual physical time. This approach significantly reduces computation time while maintaining simulation accuracy.

[0039] After iterative calculations, rutting resistance performance evaluation indices are extracted from the calculation results. The interface fusion degree index reflects the degree of diffusion and mixing between the new and old asphalt. A lower limit for mixing judgment is set to 0.2, and an upper limit is set to 0.8. The total number of grid points in the interface mixing zone is obtained by counting the number of grid points with phase field sequence parameter values ​​between 0.2 and 0.8 within the computational domain. The percentage of the total number of grid points in the interface mixing zone to the total number of grid points in the two-dimensional rectangular computational domain is used to obtain the interface fusion degree index. A higher total number of grid points in the interface mixing zone indicates a wider diffusion range and more thorough mixing between the new and old asphalt, resulting in better recycling performance. The simulated maximum rutting depth index reflects the degree of deformation of asphalt under load. Several monitoring grid points are selected vertically below the load application area, and the cumulative vertical displacement of each monitoring grid point during the iterative calculation process is recorded. The maximum cumulative vertical displacement is selected as the simulated maximum rutting depth index. A larger cumulative vertical displacement indicates more severe high-temperature flow deformation of the asphalt and poorer rutting resistance. The high-temperature flow activity index reflects the overall flow activity of asphalt under high-temperature conditions. A flow threshold of 0.005 grid units per time step is set, and the number of grid points with a velocity vector modulus greater than the flow threshold within the computational domain is counted. The percentage of this number to the total number of grid points is then calculated to obtain the high-temperature flow activity index. A higher high-temperature flow activity index indicates a large-scale active flow within the asphalt, and a greater risk of rutting.

[0040] The mix proportion parameters are adjusted based on the comparison results between the rutting resistance performance evaluation index and the preset target threshold. The target threshold for fusion is set at 15 percentage points. When the fusion index at the interface between new and old asphalt is less than 15 percentage points, it indicates insufficient diffusion mixing between the new and old asphalt. This requires increasing the warm mix additive dosage to reduce asphalt viscosity and promote diffusion, or increasing the mixing temperature to accelerate molecular thermal motion. The allowable rutting depth is set at 3 mm. When the simulated maximum rutting depth index is greater than 3 mm, it indicates insufficient high-temperature deformation resistance of the mixture. This requires reducing the milling aggregate dosage to decrease the proportion of aged asphalt, or increasing the proportion of hard components in the new asphalt to improve overall stiffness. The upper limit for flow activity is set at 8 percentage points. When the high-temperature flow activity index is greater than 8 percentage points, it indicates that the asphalt flows too actively at high temperatures. This requires increasing the dosage of polymeric modifiers in the new asphalt to enhance the elastic recovery ability of the asphalt and inhibit viscous flow. After adjusting the milling aggregate dosage, warm mix additive dosage, new asphalt dosage, and new asphalt component ratio based on the above comparison results, the process returns to step two to re-establish the calculation domain and execute simulation calculations, forming an iterative optimization loop. When the interface fusion index of new and old asphalt, the simulated maximum rutting depth index, and the high-temperature flow activity index all meet their respective target requirements, the iteration cycle is terminated, and the mix proportion parameters at this time are used as the final mix proportion parameters for the actual production of high-volume milled material recycled warm-mix asphalt mixtures.

[0041] As an optional implementation, the computational domain can also employ a three-dimensional structure to more realistically reflect the spatial characteristics of asphalt mixtures. Establishing a three-dimensional computational domain requires a D3Q19 or D3Q27 lattice structure, significantly increasing the number of lattice points and computational load, suitable for refined analysis scenarios with high accuracy requirements. In the phase field update sub-step, the Cahn-Hilliard equation can be used instead of the simple diffusion-convection equation. The Cahn-Hilliard equation naturally describes phase separation and interfacial tension effects, suitable for systems with a significant tendency for phase separation. The lattice Boltzmann collision model can also use a multi-relaxation-time model instead of a single-relaxation-time model. The multi-relaxation-time model has better numerical stability when dealing with fluids with high viscosity ratios, suitable for situations where there is a large difference in viscosity between new and old asphalt.

[0042] Taking a major highway overhaul project as an example, the specific implementation process of this invention is explained in detail. The milled material recycled in this project comes from asphalt pavement with a service life of more than 12 years. The old asphalt content in the milled material is about 4.8%. It is planned to use a large-volume milled material with 50% asphalt content to produce recycled warm-mix asphalt mixture.

[0043] First, physical property parameters were obtained. A random sample of 250 grams of milled material was taken from the milled material transported from the site and placed in the filter paper tube of the extractor. 500 ml of trichloroethylene solvent was added to the extraction flask for extraction. After approximately 3 hours of cyclic extraction, the reflux liquid became colorless and transparent. The asphalt solution in the extraction flask was transferred to a distillation apparatus, where the solvent was removed by distillation at 160 degrees Celsius, ultimately recovering approximately 12 grams of old asphalt.

[0044] Penetration tests were performed on the recycled asphalt. The asphalt was heated to 140°C until it was fluid, then poured into a standard penetration test dish. After cooling to room temperature for 1.5 hours, it was placed in a 25°C constant temperature water bath and kept warm for 1.5 hours. The test dish was then placed on a penetration tester, allowing a 100g standard needle to freely penetrate the asphalt sample within 5 seconds. Tests were conducted at three different locations on the same sample. The three test results were 28, 26, and 27, respectively. The arithmetic mean was taken as the penetration value of the recycled asphalt. The value is 27, and the unit is 0.1 millimeters.

[0045] The softening point of the used asphalt was tested. The used asphalt was poured into a standard brass sample ring and cooled to room temperature. A standard steel ball with a diameter of 9.53 mm and a mass of 3.5 g was placed in the center of the sample ring, and the entire apparatus was immersed in distilled water. The temperature was increased at a rate of 5 degrees Celsius per minute, and the temperature was recorded when the asphalt sank to contact the metal base below. Two parallel tests were conducted, with results of 62 degrees Celsius and 63 degrees Celsius respectively. The average value was taken as the softening point value of the used asphalt. It is 62.5 degrees Celsius.

[0046] Road petroleum asphalt No. 70 was selected as the new asphalt. Penetration tests were conducted on the new asphalt using the same method. The results of three tests were 68, 71, and 70, respectively. The arithmetic mean was taken as the penetration value of the new asphalt. The value is 70, and the unit is 0.1 mm. Softening point tests were conducted on new asphalt. Two parallel test results were 48 degrees Celsius and 49 degrees Celsius, respectively. The average value was taken as the softening point value of the new asphalt. The temperature is 48.5 degrees Celsius.

[0047] A surfactant-based warm mix agent with a melting point of approximately 95 degrees Celsius was selected. Eight milligrams of the warm mix agent sample were weighed and placed in an aluminum sample crucible of a differential scanning calorimeter. The sample was heated from room temperature to 150 degrees Celsius at a heating rate of 10 degrees Celsius per minute, and the heat flow curve was recorded. The melting point of the warm mix agent was obtained by reading the temperature corresponding to the onset of the endothermic peak from the heat flow curve. It is 94 degrees Celsius.

[0048] The particle size distribution of the milled material was tested using a laser particle size analyzer. The results showed that particles smaller than 0.075 mm accounted for 8.2%, particles between 0.075 mm and 2.36 mm accounted for 32.5%, particles between 2.36 mm and 13.2 mm accounted for 45.8%, and particles larger than 13.2 mm accounted for 13.5%. A 50 mm cube milled material sample was scanned using X-ray computed tomography (CT) at a resolution of 15 micrometers. The initial porosity of the milled material was extracted from the acquired three-dimensional spatial distribution image, yielding a value of 6.3% and an average pore diameter of 0.8 mm.

[0049] The simulation then proceeds to the coupled simulation stage using the phase-field method and the lattice Boltzmann method. A two-dimensional rectangular computational domain is established in the computer memory, with 500 grid points horizontally and 200 grid points vertically. Each grid point represents a physical size of 1 micrometer. The computational domain is divided into three segments horizontally: the first segment occupies the area from the 1st to the 150th grid point (150 grid points); the second segment occupies the area from the 151st to the 350th grid point (200 grid points); and the third segment occupies the area from the 351st to the 500th grid point (150 grid points).

[0050] Define phase-field sequence parameters Characterizes the asphalt phase properties at any location within the computational domain. This applies to all grid points within the first computational segment. Assign an initial value of 0 to all grid points within the third calculation segment. Assign an initial value of 1. For grid points within the second calculation segment, located in the horizontal direction... At each grid point Initial value according to Calculation, where This is the horizontal index number of the grid point. For example, at grid point 151... At the 250th grid point At the 350th grid point .

[0051] The reference viscosity value of old asphalt is calculated based on its penetration and softening point values. An empirical correlation formula is used. Perform calculations, where Dynamic viscosity, measured in Pascal-seconds. Temperature, in degrees Celsius. This is an empirical coefficient related to penetration and softening point. For old asphalt with a penetration value of 27 and a softening point of 62.5 degrees Celsius, the reference viscosity value of the old asphalt at a reference temperature of 60 degrees Celsius is... The calculated value is approximately 850 Pascals per second. For new asphalt with a penetration value of 70 and a softening point of 48.5 degrees Celsius, the reference viscosity value for new asphalt at a reference temperature of 60 degrees Celsius is... The calculation yields approximately 180 Pascals per second.

[0052] Set the simulated temperature value The temperature was 135 degrees Celsius, which is higher than the softening point of old asphalt (62.5 degrees Celsius), higher than the softening point of new asphalt (48.5 degrees Celsius), and higher than the melting point of warm mix additives (94 degrees Celsius). Based on the simulated temperature value, the reference viscosity value was corrected for temperature. The asphalt viscosity decreased exponentially with increasing temperature. The corrected viscosity value for old asphalt was... Approximately 45 Pascals per second, the new asphalt viscosity value after temperature correction. It takes approximately 6 Pascals of a second.

[0053] A lattice Boltzmann velocity discrete lattice is constructed, using a D2Q9 lattice structure. Nine discrete velocity vectors are defined, with the velocity vector in the 0th direction being... The first direction velocity vector The second direction velocity vector The third direction velocity vector The fourth direction velocity vector The fifth direction velocity vector The 6th direction velocity vector The 7th direction velocity vector The 8th direction velocity vector Nine particle distribution functions are initialized at each grid point. to The initial values ​​are all set to the equilibrium distribution function values ​​in the corresponding directions, and the initial macroscopic density is... Set to 1, initial flow velocity vector is set to .

[0054] Establish the correlation between phase-field order parameters and local kinematic viscosity. For any grid point, read the phase-field order parameter value at the current grid point. Local kinematic viscosity at the current grid point Calculated using linear interpolation. For example, when hour Pascal second, when hour Pascal second, when hour Pascal-second. The dynamic viscosity is converted to grid units. The conversion factor is determined based on the principle of similarity. One grid unit of viscosity corresponds to an actual viscosity of 10 Pascal-seconds. Therefore, the viscosity of old asphalt in grid units is 4.5, and the viscosity of new asphalt is 0.6.

[0055] Set boundary conditions. The upper and lower boundaries use periodic boundary conditions: when the particle distribution function migrates out from row 200 of the upper boundary, it migrates in from row 1 of the lower boundary; conversely, when the particle distribution function migrates out from row 1 of the lower boundary, it migrates in from row 200 of the upper boundary. The left boundary uses a constant velocity inlet boundary condition, with the inlet velocity value... The time step is set to 0.01 grid units. The right boundary uses a constant pressure outlet boundary condition, and the density value corresponding to the outlet pressure is set to 1.

[0056] External load conditions simulating rut formation are applied. The load application region is selected at the 200th to 300th grid points horizontally along the upper boundary, with a width of 100 grid points. The simulated wheel load pressure is set to 0.7 MPa, which translates to a volumetric force source term of -0.0001 grid units (the negative sign indicates a downward direction). The particle distribution function along the vertically negative direction is calculated at each grid point within the load application region. Superimposed volumetric force source terms.

[0057] Perform alternating iterative calculations, setting the total number of iterations to 40,000. The phase field update sub-step is illustrated using the calculation process at the 250th grid point horizontally and the 100th grid point vertically at step 20,000 as an example. Read the current phase field sequence parameter value of the target grid point. Read the phase field sequence parameter value of the adjacent grid point in the positive horizontal direction, i.e., the 251st grid point. Read the phase field sequence parameter value of the adjacent grid point in the negative horizontal direction, i.e., the 249th grid point. Read the phase field sequence parameter values ​​of the adjacent grid points in the vertical positive direction, i.e., the 101st row. Read the phase field sequence parameter value of the adjacent grid point in the vertical negative direction, i.e., the 99th row. The arithmetic mean of the four adjacent phase field sequence parameter values ​​is used to obtain the adjacent average phase field sequence parameter value. Calculate the phase field diffusion driving force. .

[0058] Obtain the current flow velocity vector at the target grid point, including the horizontal component. Vertical component Calculate horizontal convective transport volume. Calculate vertical convective transport volume. Calculate the total convective transport volume. Set the phase-field diffusion coefficient. The diffusion update increment is calculated as follows: Setting the convective transport coefficient to 0.5, the convective update increment is calculated as follows: Calculate the updated phase-field sequence parameter values. This value is between 0 and 1 and does not require correction.

[0059] The calculation process of the lattice Boltzmann collision migration sub-step will be further explained. For the same target lattice point, the current macroscopic density value is first calculated. The sum of the nine particle distribution function values ​​is obtained. Calculate the components of the current flow velocity vector. , ,in and The first The horizontal and vertical components of each discrete velocity direction are calculated. , .

[0060] Perform a collision operation on each particle distribution function. Taking the first direction, i.e., the positive horizontal direction, as an example, calculate the local equilibrium distribution function value. Weighting coefficients for each direction in the D2Q9 lattice structure. for: , , The formula for calculating the local equilibrium distribution function is: ,in For the first Discrete velocity vectors in each direction, This is the macroscopic velocity vector. For the first direction, , Substituting into the calculation, we get... Read the current particle distribution function value. Calculate the deviation from equilibrium. Based on the phase-field sequence parameter values ​​at the target grid point. Calculate local kinematic viscosity Grid unit. Relaxation time value calculated based on local kinematic viscosity. The collision adjustment amount is calculated as follows: Calculate the particle distribution function value after the collision. .

[0061] The collision calculation process is repeated for the particle distribution functions in the remaining eight directions. After the collision operation is completed, a migration operation is performed to move the target grid point after the collision. Move to the adjacent grid point in the positive horizontal direction, i.e., the 251st grid point. Move to the adjacent grid point in the vertical positive direction, i.e., row 101. Move to the adjacent grid point in the negative horizontal direction, i.e., the 249th grid point. Move to the adjacent grid point in the negative vertical direction, i.e., row 99. Move to the adjacent grid point in the upper right corner, and Move to the adjacent grid point in the upper left corner, and Move to the adjacent grid point in the lower left corner, and Move to the adjacent grid point in the lower right corner. The particles are retained at the target grid point. After migration, the received particle distribution function is summarized at each grid point, and the macroscopic density and velocity vector are recalculated.

[0062] After a total of 40,000 iterations, the rutting resistance performance evaluation index was extracted. A lower limit for mixing was set to 0.2, and an upper limit to 0.8. All 100,000 grid points in the computational domain were traversed, and 11,850 grid points had phase field order parameter values ​​between 0.2 and 0.8. The fusion index of the new and old asphalt interface was calculated. 101 monitoring grid points were selected directly below the load application area, located at rows 50 to 150 vertically and the 250th horizontally. The cumulative vertical displacement of each grid point was recorded. The maximum value of 3.72 mm was observed at row 80, which is the simulated maximum rut depth index. A flow threshold of 0.005 grid units per time step was set. The number of grid points with a flow velocity vector magnitude greater than 0.005 was 9230. The high-temperature flow activity index was calculated as follows: .

[0063] Entering the mix proportion parameter adjustment stage. The target threshold for fusion degree is set at 15 percentage points. Currently, the fusion degree index at the interface between new and old asphalt is 11.85 percentage points, which is less than 15 percentage points, indicating insufficient fusion at the interface. The allowable rutting depth is set at 3 mm. Currently, the simulated maximum rutting depth index is 3.72 mm, which is greater than 3 mm, indicating insufficient rutting resistance. The upper limit for flow activity is set at 8 percentage points. Currently, the high-temperature flow activity index is 9.23 percentage points, which is greater than 8 percentage points, indicating excessively active high-temperature flow.

[0064] Based on the above assessment results, the mix proportions were adjusted. To address the issue of insufficient interfacial fusion, the warm mix admixture dosage was increased from 3% to 4% of the asphalt mass. To address the issue of insufficient rutting resistance, the milling aggregate dosage was reduced from 50% to 45%. To address the issue of excessively active high-temperature flow, a styrene-butadiene-styrene block copolymer modifier was added to the new asphalt at a dosage of 4% of the new asphalt mass.

[0065] Returning to step two, the simulation calculation was re-executed using the adjusted mix proportions. Due to the addition of the modifier, the viscosity of the new asphalt increased. The penetration value of the modified new asphalt was re-measured to be 58, and the softening point to be 56 degrees Celsius. The calculated temperature-corrected viscosity of the new asphalt was approximately 12 Pascals per second. After updating the correlation parameters between the phase field sequence parameters and the local kinematic viscosity, the 40,000-step iterative calculation was re-executed.

[0066] After the second round of simulations, evaluation indicators for rutting resistance were extracted. The number of grid points with phase field order parameter values ​​between 0.2 and 0.8 was 16,320, and the calculated fusion degree index of the new and old asphalt interface was 16.32 percentage points, exceeding the target fusion degree threshold by 15 percentage points. The maximum cumulative vertical displacement of the monitored grid points was 2.65 mm, which corresponds to a simulated maximum rutting depth of 2.65 mm, less than the allowable rutting depth of 3 mm. The number of grid points with flow velocity vector modulus values ​​greater than 0.005 was 6,850, and the calculated high-temperature flow activity index was 6.85 percentage points, less than the upper limit of flow activity by 8 percentage points.

[0067] When the interface blending degree of new and old asphalt, the simulated maximum rutting depth, and the high-temperature flow activity all meet their respective target requirements, the iterative optimization cycle is terminated. The mix proportions at this point are then determined as the final mix proportions: milling aggregate content is 45 percentage points, warm mix additive content is 4 percentage points of the asphalt mass, and the new asphalt is modified asphalt with 4 percentage points of styrene-butadiene-styrene block copolymer added; the new asphalt content is determined based on the asphalt-aggregate ratio requirements. Actual production of high-volume milling aggregate recycled warm mix asphalt mixtures is conducted using these final mix proportions, with the production temperature controlled between 135°C and 145°C.

[0068] refer to Figure 3 The horizontal axis of the graph represents the number of iterations in the LBM simulation. The unit is time step, which has a linear mapping relationship with actual physical time and is used to characterize the time dimension of the simulation process. The left vertical axis represents the simulated maximum rut depth. The unit is millimeters; the right vertical axis represents the degree of integration between the old and new asphalt interfaces. The units are percentages. The figure contains three curves that evolve over time, corresponding to the performance response and interface behavior under different mix proportions. The first curve, the red dashed line corresponding to the left vertical axis, represents the evolution trajectory of rut depth under initial mix proportions (e.g., high milling mix content or insufficient warm mix agent). With each iteration... As the load increases, the curve shows a rapid upward trend, indicating that the asphalt mixture has undergone significant cumulative permanent deformation under continuous load, and no obvious signs of stable convergence have appeared in the later stage of the simulation, showing poor rutting resistance.

[0069] The permissible threshold for rut depth is marked by a horizontal red dashed line in the diagram. For example, at a depth of 3 mm, the initial mix design curve clearly exceeded this threshold. The second curve, the solid red line corresponding to the left ordinate, represents the evolution trajectory of rut depth after mix design optimization (e.g., adjusting the ratio of new and old asphalt or adding modifiers). The upward slope of this curve is significantly less than that of the initial mix design curve, and it gradually flattens out in the later stages of the simulation, eventually converging to... The following figures demonstrate that the optimized asphalt mixture exhibits higher stiffness and better elastic recovery, effectively resisting high-temperature deformation. The third curve, corresponding to the blue dotted line on the right-hand vertical axis, represents the degree of integration between the new and old asphalt interfaces. The curve exhibits an evolutionary pattern over time. It roughly follows the square root time law or a similar diffusion growth pattern, increasing with the number of iterations. The increase, The value gradually increases, reflecting the continuous diffusion and fusion of new and old asphalt molecules under the action of thermal motion and mechanical stirring. The target fusion threshold is marked by a horizontal blue dashed line in the figure. For example, 15%. The curve is related to... The number of iterations corresponding to the intersection points indicates the shortest mixing time required to achieve the desired regeneration effect. Figure 3 By comparing mechanical performance indicators (rutting depth) and physicochemical indicators (interfacial fusion) on the same time axis, the decision-making logic of the control method is clearly illustrated: that is, to find an optimal combination of ratio parameters that ensures the interfacial fusion. Exceed At the same time, the maximum rut depth Controlled Within this range, the dynamic evolution curve provides a quantitative analytical tool for determining the final production process parameters.

[0070] The present invention has been described in detail above. Specific examples have been used to illustrate the principles and implementation methods of the invention. The descriptions of the embodiments above are merely for the purpose of helping to understand the method and core ideas of the present invention. It should be noted that those skilled in the art can make various improvements and modifications to the present invention without departing from its principles, and these improvements and modifications also fall within the protection scope of the claims of the present invention.

Claims

1. A method for controlling rutting resistance in warm-mix asphalt with high admixture of milled aggregate, characterized in that, Includes the following steps: Step 1: Conduct penetration and softening point tests on the old asphalt in the milled material to obtain the penetration value and softening point value of the old asphalt. Conduct penetration and softening point tests on the new asphalt to be added to obtain the penetration value and softening point value of the new asphalt. Conduct melting point tests on the warm mix additive to be added to obtain the melting point value of the warm mix additive. Step 2: Establish a two-dimensional rectangular computational domain and divide it horizontally into a first computational segment, a second computational segment, and a third computational segment. The first computational segment represents the old asphalt phase region, the third computational segment represents the new asphalt phase region, and the second computational segment represents the transition region between the old and new asphalt interfaces. Establish the phase field sequence parameter distribution to represent the phase state properties of asphalt, and establish a lattice Boltzmann velocity discrete lattice to simulate asphalt flow behavior. Perform alternating iterative calculations of the phase field update sub-step and the lattice Boltzmann collision migration sub-step, apply external load conditions simulating rutting formation, and extract anti-rutting performance evaluation indicators after iteration. Step 3: Compare the anti-rutting performance evaluation index with the preset target threshold, adjust the mix proportion parameters according to the comparison results, and return to Step 2 to re-simulate until the target requirements are met, and determine the final mix proportion parameters.

2. The method according to claim 1, characterized in that, Step one also includes: using a laser particle size analyzer to test the particle size distribution of aggregates in the milled material to obtain aggregate particle size distribution data; using an X-ray computed tomography (CT) scanner to scan the milled material sample to obtain a three-dimensional spatial distribution image of the internal pores of the milled material, and extracting the initial porosity value and the average pore diameter value of the milled material from the three-dimensional spatial distribution image.

3. The method according to claim 1, characterized in that, In step two, the phase field sequence parameter ranges from zero to one in a continuous interval. When the phase field sequence parameter is zero, it indicates that the corresponding position is entirely old asphalt phase. When the phase field sequence parameter is one, it indicates that the corresponding position is entirely new asphalt phase. When the phase field sequence parameter is between zero and one, it indicates that the corresponding position is in a mixed transition state of old and new asphalt. When initializing the phase field sequence parameter, the initial value of the phase field sequence parameter at all positions in the first calculation section is set to zero, the initial value of the phase field sequence parameter at all positions in the third calculation section is set to one, and the initial value of the phase field sequence parameter at each position in the second calculation section is set to a transitional distribution that gradually changes from zero to one in the horizontal direction.

4. The method according to claim 1, characterized in that, In step two, the Boltzmann velocity discretization lattice uses a D2Q9 lattice structure for velocity space discretization. The D2Q9 lattice structure sets nine discrete velocity directions at each lattice point. The nine discrete velocity directions include one zero velocity direction, four unit velocity directions along the positive and negative coordinate axes, and four unit velocity directions along the diagonal directions. Nine particle distribution functions are defined at each lattice point, and the nine particle distribution functions correspond to the nine discrete velocity directions. Each particle distribution function characterizes the number density of virtual particles moving along the corresponding discrete velocity direction.

5. The method according to claim 1, characterized in that, Step two also includes: setting a simulated temperature value that is higher than the softening point of the old asphalt, higher than the softening point of the new asphalt, and higher than the melting point of the warm mix agent; establishing a correlation rule between the phase field sequence parameter and the lattice Boltzmann local viscosity; determining the reference viscosity value of the old asphalt based on its penetration and softening point; determining the reference viscosity value of the new asphalt based on its penetration and softening point; and performing temperature correction on the reference viscosity values ​​of the old and new asphalt based on the simulated temperature value to obtain the temperature-corrected viscosity values ​​of the old and new asphalt. For any grid point within the computational domain, when the current phase field sequence parameter at the current grid point is zero, the local kinematic viscosity at the current grid point is set to the temperature-corrected old asphalt viscosity value. When the current phase field sequence parameter at the current grid point is one, the local kinematic viscosity at the current grid point is set to the temperature-corrected new asphalt viscosity value. When the current phase field sequence parameter at the current grid point is between zero and one, the local kinematic viscosity at the current grid point is set to the linear interpolation value between the temperature-corrected old asphalt viscosity value and the temperature-corrected new asphalt viscosity value, with the interpolation coefficients of the linear interpolation being the current phase field sequence parameter value.

6. The method according to claim 1, characterized in that, The execution process of the phase field update sub-step in step two is as follows: For any target grid point in the computational domain, obtain the current phase field sequence parameter value of the target grid point, obtain the current phase field sequence parameter values ​​of the target grid point in the horizontal positive direction, the horizontal negative direction, the vertical positive direction, and the vertical negative direction, a total of four adjacent phase field sequence parameter values, calculate the arithmetic mean of the four adjacent phase field sequence parameter values ​​to obtain the adjacent average phase field sequence parameter value, and calculate the difference between the adjacent average phase field sequence parameter value and the current phase field sequence parameter value of the target grid point to obtain the phase field diffusion driving amount; Obtain the current velocity vector at the target grid point. Calculate the horizontal convective transport by multiplying the horizontal component of the current velocity vector by the difference between the current phase sequence parameter values ​​of adjacent grid points in the positive and negative horizontal directions. Calculate the vertical convective transport by multiplying the vertical component of the current velocity vector by the difference between the current phase sequence parameter values ​​of adjacent grid points in the positive and negative vertical directions. Pair the horizontal and vertical convective transport values. The total convective transport is obtained by adding the phase field diffusion driving quantity and multiplying it by the preset phase field diffusion coefficient to obtain the diffusion update increment. The total convective transport is multiplied by the preset convective transport coefficient to obtain the convective update increment. The current phase field sequence parameter value of the target grid point is added to the diffusion update increment and then subtracted from the convective update increment to obtain the updated phase field sequence parameter value. When the updated phase field sequence parameter value is less than zero, the updated phase field sequence parameter value is corrected to zero. When the updated phase field sequence parameter value is greater than one, the updated phase field sequence parameter value is corrected to one.

7. The method according to claim 4, characterized in that, The execution process of the lattice Boltzmann collision migration sub-step in step two is as follows: For any target lattice point within the computational domain, a collision operation is first performed. For each particle distribution function at the target lattice point, the local equilibrium distribution function value corresponding to the current particle distribution function is calculated based on the current macroscopic density value and the current velocity vector at the target lattice point. The difference between the current particle distribution function value and the local equilibrium distribution function value is calculated to obtain the deviation from equilibrium. The relaxation time value is calculated based on the local kinematic viscosity value at the target lattice point. The deviation from equilibrium is divided by the relaxation time value to obtain the collision adjustment amount. The current particle distribution function value is then subtracted from the value of the local kinematic viscosity value. The collision adjustment amount yields the particle distribution function value after the collision. After the collision operation is completed, a migration operation is performed, migrating the nine particle distribution functions after the collision at the target grid point to adjacent grid points along their respective discrete velocity directions. The particle distribution function after the collision corresponding to the zero velocity direction is retained at the target grid point. After the migration operation is completed, the particle distribution functions migrated to the current grid point are summarized for each grid point. The nine particle distribution function values ​​are added together to obtain the updated macroscopic density value. The nine particle distribution function values ​​are multiplied by their respective discrete velocity vectors, added together, and then divided by the updated macroscopic density value to obtain the updated flow velocity vector.

8. The method according to claim 1, characterized in that, Step two also includes setting boundary conditions at the boundaries of the two-dimensional rectangular computational domain: setting periodic boundary conditions at the upper and lower boundaries of the two-dimensional rectangular computational domain such that the particle distribution function migrating out from the upper boundary migrates into the lower boundary and the particle distribution function migrating out from the lower boundary migrates into the upper boundary; setting constant flow velocity inlet boundary conditions at the left boundary of the two-dimensional rectangular computational domain; and setting constant pressure outlet boundary conditions at the right boundary of the two-dimensional rectangular computational domain.

9. The method according to claim 1, characterized in that, The method for applying the simulated rutting external load conditions in step two is as follows: the central section of the upper boundary of the two-dimensional rectangular computational domain is selected as the load application area, and a downward volumetric force source term is superimposed on the particle distribution function at each grid point within the load application area; the rutting resistance performance evaluation index includes the new and old asphalt interface fusion index, the simulated maximum rutting depth index, and the high-temperature flow activity index. The new and old asphalt interface fusion index is the percentage of grid points whose phase field sequence parameter values ​​are between the preset lower limit and upper limit of the mixing judgment value to the total number of grid points in the two-dimensional rectangular computational domain. The simulated maximum rutting depth index is the maximum value of the cumulative vertical displacement of the monitoring grid points directly below the load application area. The high-temperature flow activity index is the percentage of grid points whose velocity vector modulus is greater than the preset flow threshold to the total number of grid points in the two-dimensional rectangular computational domain.

10. The method according to claim 9, characterized in that, The specific methods for adjusting the proportioning parameters in step three are as follows: when the fusion index of the new and old asphalt interface is less than the preset fusion target threshold, increase the amount of warm mix additive or increase the mixing temperature; when the simulated maximum rutting depth index is greater than the preset allowable rutting depth value, reduce the amount of milling material or increase the proportion of hard components in the new asphalt; when the high temperature flow activity index is greater than the preset upper limit of flow activity, increase the amount of polymer modifier in the new asphalt.

Images