A Discrete Element Method for Evaluating the Contact Anisotropy of Asphalt Mixture Skeleton

Abstract

This invention discloses a discrete element method for evaluating the contact anisotropy of asphalt mixture skeletons, comprising the following steps: (1) selecting limestone aggregates of different particle size ranges and performing CT scanning and 3D modeling reconstruction; (2) determining the constitutive model parameters between aggregates and asphalt mortar through uniaxial creep tests, and generating virtual samples of asphalt mixtures according to the gradation type; (3) classifying the contact normals and normal contact forces of the above contacts according to the angle between the contact vector and the coordinate axis; (4) converting the above contact normal vectors into unit vectors, calculating the fabric tensor of the contact normals, and converting the contact normals and normal contact force vectors into polar coordinates, and drawing wind rose diagrams at specified angular intervals; (5) constructing Fourier series functions and fitting the wind rose diagrams of the contact normals and normal contact forces, and evaluating the anisotropy of asphalt mixture skeleton contact with the fitted anisotropy parameters and fabric tensor.

Description

Technical Field

[0001] This invention relates to discrete element numerical simulation and skeleton evaluation technology for asphalt mixtures, and particularly to a method for evaluating the contact anisotropy of asphalt mixture skeletons based on discrete element analysis. Background Technology

[0002] Asphalt mixtures are multiphase heterogeneous composite materials containing coarse aggregates, fine aggregates, mineral fillers, asphalt binders, and voids. Their mechanical properties are directly related to the spatial structure of each component, such as the distribution of coarse aggregates, the skeleton formed by the aggregate particles, and the size of the voids. The skeleton, composed of several adjacent aggregates, has the ability to resist and transfer external load stresses. By analyzing the dynamic response of the aggregate skeleton, a comprehensive understanding of the structural characteristics and service performance of asphalt mixtures can be obtained.

[0003] Previously, anisotropy of asphalt mixtures could be measured by examining modulus differences in different directions. Furthermore, digital image processing techniques could be used to characterize anisotropy through morphological features such as the spatial distribution of aggregates and voids within the sample and the orientation of the long axis. Compared to experimental methods, increasing research focuses on numerical methods to simulate the skeleton contact behavior of asphalt mixtures during deformation. In particular, the Discrete Element Method (DEM) is more suitable than the continuum-based finite element method for analyzing the microscopic skeleton contact behavior of mixed particles. By using Particle Flow Code (PFC) and other custom modeling platforms, simulations of real aggregates can be achieved, capturing the motion of individual aggregates. Virtual testing on asphalt mixture models composed of these aggregates provides a promising approach for studying contact force chains and contact anisotropy. Aggregate properties, temperature, and mixture type significantly influence the contact forces and anisotropy of the skeleton structure under external loads, with asphalt mortar and aggregate skeleton playing crucial roles in the rutting resistance of asphalt mixtures.

[0004] Undeniably, significant progress has been made in using DEM (Depth-Earth Model) to describe the skeleton contact behavior between coarse aggregates in asphalt mixtures and to elucidate the rutting-related mechanisms. The interactions between aggregates and their contact anisotropy under load are crucial for evaluating the rutting resistance and long-term service behavior of asphalt mixtures. Extracting key features through DEM analysis allows for a deeper understanding of the complex relationship between specific structural configurations and pavement durability. Determining the spatial orientation and contact behavior of coarse aggregates and asphalt mortar can reveal the intrinsic mechanisms of skeleton structure evolution under load. Therefore, it is necessary to propose a method for evaluating the skeleton contact anisotropy of asphalt mixtures to establish the relationship between contact anisotropy and mechanical properties. Summary of the Invention

[0005] Objective of the Invention: To address the above-mentioned problems, the objective of this invention is to provide a discrete element method for evaluating the contact anisotropy of asphalt mixture skeletons, effectively assessing the contact anisotropy of asphalt mixtures under load. This method helps overcome the limitations of previous methods that solely evaluated the degree of anisotropy based on differences in modulus in different directions and aggregate distribution. By understanding the spatial orientation of the aggregate-asphalt mortar contact during rutting, a theoretical basis can be provided for establishing the relationship between contact anisotropy and mechanical properties, and clarifying the rutting development mechanism of asphalt mixtures.

[0006] Technical solution: The present invention provides a method for evaluating the contact anisotropy of asphalt mixture skeleton based on discrete element method, comprising the following steps:

[0007] (1) Select limestone aggregates with different particle size ranges and perform CT scanning and 3D modeling reconstruction;

[0008] (2) Import the above aggregate coordinate file into the discrete element software, determine the constitutive model parameters between aggregate and asphalt mortar through uniaxial creep test, and generate virtual samples of asphalt mixture according to the gradation type;

[0009] (3) Apply load to virtual sample, identify contact between aggregate and asphalt mortar, and classify the contact normal and normal contact force of the above contact according to the angle between the contact vector and the coordinate axis.

[0010] (4) Convert the above contact normal vector into a unit vector, calculate the structure tensor of the contact normal, and convert the contact normal and normal contact force vector into polar coordinates, and draw the wind rose diagram at specified angle intervals.

[0011] (5) Construct Fourier series functions and fit wind rose diagrams of contact normals and normal contact forces to evaluate the anisotropy of asphalt mixture skeleton contact using the fitted anisotropy parameters and fabric tensors.

[0012] Furthermore, in step (1), selecting limestone aggregates with different particle size ranges and performing CT scanning and 3D modeling reconstruction includes:

[0013] (101) The aggregate obtained from the quarry is screened and classified into four categories according to the particle size based on the screen aperture size: 2.36-4.75mm, 4.75-9.5mm, 9.5-13.2mm and 13.2-16.0mm.

[0014] (102) At least 25 aggregates of each particle size were selected, washed and fixed in a box. An industrial computer tomography scanner was used to scan the aggregates in the box from front to back to obtain the tomographic projection image of the aggregates.

[0015] (103) The scanned tomographic images were imported into AVIZO software for morphological processing, including filtering, noise reduction, contrast enhancement and other methods. Threshold segmentation was used for binarization to identify the planar regions containing aggregates. The three-dimensional structural models of each aggregate were obtained by using the software's automatic stacking.

[0016] (104) Continue to generate a surface mesh covering the three-dimensional aggregate profile in the software, and simplify the surface mesh of the aggregate by controlling the number of points, surfaces and edges.

[0017] Furthermore, in step (2), the above aggregate coordinate file is imported into the discrete element method software, and the constitutive model parameters between the aggregate and the asphalt mortar are determined through uniaxial creep tests. Virtual samples of the asphalt mixture are then generated according to the gradation type, including:

[0018] (201) After simplifying the surface mesh of the aggregate, its contour coordinate information is output to an STL file, and a command is written to import the above STL file into the Discrete Element Method (PFC) software to create an aggregate cluster template. By adjusting the two parameters that affect the simulation accuracy of the aggregate cluster template, namely distance and ratio, a balance between simulation accuracy and simulation efficiency is achieved;

[0019] (202) The elastic modulus of the aggregates was obtained by nanoindentation test, and the constitutive model parameters of the contact between the aggregates were determined by selecting a linear stiffness model. The expression is as follows:

[0020] k n =2EL

[0021]

[0022] In the formula, E is the elastic modulus of the aggregate, which is determined to be 55.5 MPa according to the nanoindentation test; L is the distance between the two contact ball elements, L=R1+R2, which is 2 mm in this study; υ is the Poisson's ratio of the aggregate, which is 0.25;

[0023] (203) After constructing the aggregate model, asphalt mortar containing asphalt binder and fine aggregate with a particle size of less than 2.36 mm was prepared in the laboratory and uniaxial creep loading tests were conducted at different temperatures. The creep model parameters were obtained by fitting the test data using Burger's viscoelastic model, and its expression is as follows:

[0024]

[0025] In the formula, E1 is the elastic modulus of the Kelvin portion of the Burger's model; η1 is the viscosity of the Kelvin portion of the Burger's model; E2 is the elastic modulus of the Maxwell portion of the Burger's model; η2 is the viscosity of the Maxwell portion of the Burger's model; and t is the loading time.

[0026] Simultaneously, the creep model parameters obtained above are converted into constitutive model parameters for the contact between aggregates and asphalt mortar in discrete element software, as expressed below:

[0027] C mn =η1L

[0028] K mn =E1L

[0029] C kn =η2L

[0030] K kn =E2L

[0031]

[0032]

[0033]

[0034]

[0035] In the formula, C mn K mn C kn and K kn To simulate the parameters of the Burger's model at the normal direction; C ms K ms C ks and K ks υ′ represents the parameters of the Burger's model in the shear direction at the simulation scale; υ′ represents the Poisson's ratio of the asphalt mortar.

[0036] (204) Generate a cuboid box and corresponding six walls based on the size of the rut sample. Calculate the number of four aggregate particles in step (101) according to the required gradation type and put them into the virtual box formed by the rut sample. Scaling and expanding these aggregates in sequence to eliminate the stress caused by particle overlap.

[0037] (205) Generate a specified number of small balls with a diameter of 1.0 mm in the remaining space of the rut sample box to represent asphalt mortar, and delete a fixed number of small balls according to the required porosity of the sample to finally form a virtual rut sample of asphalt mixture.

[0038] Furthermore, in step (3), a load is applied to the virtual sample to identify the contact between the aggregate and the asphalt mortar. The contact normal and normal contact force of the above contact are classified according to the angle between the contact vector and the coordinate axis, including:

[0039] (301) After creating the virtual sample, remove the wall above the sample while keeping the walls on the sides and bottom, and generate a loading plate with a size of 50mm×50mm directly above the sample.

[0040] (302) Assign the contact constitutive model parameters to each component, and convert the moving loading time specified in the laboratory rutting test into static loading and determine the equivalent loading time, as shown in the following expression:

[0041]

[0042]

[0043] In the formula, l is the total contact length between the rubber wheel and the sample; m is the mass of the rubber wheel, which is 78 kg; g = 9.8 m / s² 2 , where is the acceleration due to gravity; p is the contact pressure, equal to 0.7 MPa; d = 0.05 m, which is the width of the rubber wheel; t is the cumulative loading time.

[0044] The constitutive model parameters obtained in step (203) were scaled by a factor of 10,000 using the time-temperature equivalence principle, and a program was written to input them into the discrete element PFC software;

[0045] (303) Use the FISH language in PFC software to activate the servo system to apply the load, and record the particle contact between aggregates, between asphalt mortars and between aggregates and asphalt mortars during loading;

[0046] (304) After loading is complete, output particle contact information, including normal contact force f. n and the three-dimensional contact normal vector n i Simultaneously, the three-dimensional contact vector is projected onto the XOZ plane to obtain the two-dimensional projection vector n of the contact normal vector. ip ;

[0047] (305) Calculate the angle θ between each contact vector and the X-axis, ranging from 0° to 360°, and sort these vectors from low to high according to the size of the angle.

[0048] Further, in step (4), the contact normal vector is converted into a unit vector, the configuration tensor of the contact normal is calculated, and the contact normal and normal contact force vectors are converted into polar coordinates. The wind rose diagram is drawn at specified angular intervals, including:

[0049] (401) Convert the three-dimensional contact normal vector obtained in step (304) into a unit vector, and calculate the configuration tensor of the contact normal, which is expressed as:

[0050]

[0051] In the formula, is the structural tensor; N is the total number of contact normal vectors in the asphalt mixture; is the component of the unit vector of the k-th contact normal; It is the tensor product;

[0052] Taking the contact normal as an example, we define the contact normal probability density function E(n), where the average value of E(n) across all directions in the entire spatial domain is 1, i.e.

[0053] ∫ Ω E(n)dΩ=1

[0054] In the formula, Ω is the solid angle corresponding to the entire surface of a unit sphere.

[0055] Therefore, the discrete contact in the above structural tensor calculation equation can be replaced by certain continuous distributions, expressed as:

[0056]

[0057] E(n) can be expanded using the second-order Fourier method, as follows:

[0058]

[0059] In the formula, D ij yes The deviation portion is represented as:

[0060]

[0061] In the formula, δ ij It's the Kronecker function:

[0062]

[0063] (402) The result is a 3×3 symmetric matrix. The three eigenvalues ​​of the construct tensor can be further extracted ( and The degree of anisotropy is measured by a tensor called the generalized octahedral structure tensor Φ.

[0064]

[0065] (403) The two-dimensional projection vector n obtained in step (304) is... ip The Cartesian coordinates were converted to polar coordinates, and wind rose diagrams were drawn at 10° angular intervals. In the rose diagram of the contact normal, the number of each bin was normalized by the total number of contacts, while in the rose diagram of the normal contact force, the number of each bin was normalized by the average normal contact force of all contacts.

[0066] Furthermore, in step (5), a Fourier series function is constructed and a wind rose diagram of the contact normal and normal contact force is fitted. The anisotropy parameters and fabric tensor obtained from the fitting are used to evaluate the anisotropy of the asphalt mixture skeleton contact, including:

[0067] (501) For the two-dimensional wind rose diagram drawn in step (403), the contact normal density function has a two-dimensional similarity in polar coordinates, namely

[0068]

[0069] In the formula, a ij It is the coefficient of contact normal anisotropy, equal to:

[0070]

[0071] Distribution function E(θ) and These are statistics of the contact normal and mean normal contact force in polar coordinates. Typically, they can be approximated by a second-order Fourier series:

[0072]

[0073]

[0074] In the formula, a is the second-order anisotropy parameter of the contact normal; θ a The principal direction that touches the anisotropy of the normal; a is the average normal contact force of all contacts in the asphalt mixture; n θ is the second-order anisotropy parameter of the normal contact force; f The principal direction of the anisotropy of the normal contact force;

[0075] (502) In the engineering drawing software OriginPro, the distribution function E(θ) and The two-dimensional wind rose diagram drawn in the fitting (403) step can be obtained by fitting a and a. n Evaluate the two-dimensional anisotropy of the contact normal and normal contact forces;

[0076] (503) The three-dimensional anisotropy of the asphalt mixture skeleton contact is evaluated using the generalized octahedral contact normal tensor Φ obtained in step (402).

[0077] Beneficial effects: Compared with the prior art, the significant advantages of this invention are:

[0078] 1. The skeleton contact anisotropy evaluation method provided by the present invention can effectively evaluate the contact anisotropy of asphalt mixture under load, and overcomes the limitations of the previous method of simply evaluating the degree of anisotropy based on the modulus difference and aggregate distribution of the sample in different directions.

[0079] 2. The discrete element method-based numerical simulation can fully explore the influencing factors of contact anisotropy in asphalt mixtures without the need for cumbersome experimental procedures.

[0080] 3. By understanding the spatial orientation of aggregates and asphalt mortar during rutting, we can provide a theoretical basis for establishing the relationship between contact anisotropy and mechanical properties and clarifying the rutting development mechanism of asphalt mixtures. Attached Figure Description

[0081] Figure 1 This is a flowchart of the present invention;

[0082] Figure 2 This is a schematic diagram of the material CT image acquisition and surface mesh simplification process in this invention;

[0083] Figure 3 This is a schematic diagram of different distance and ratio parameters in the aggregate cluster generation process of this invention;

[0084] Figure 4 This is a schematic diagram of a virtual rut sample of asphalt mixture in this invention;

[0085] Figure 5 This is a schematic diagram of the asphalt mixture loading plate in this invention;

[0086] Figure 6 This is a schematic diagram of the contact normal of the asphalt mixture and the Fourier series fitting in this invention;

[0087] Figure 7 This is a schematic diagram of the normal contact force of asphalt mixture and Fourier series fitting in this invention. Detailed Implementation

[0088] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0089] This embodiment describes a discrete element method for evaluating the contact anisotropy of asphalt mixture skeletons, comprising the following steps, the process of which is as follows: Figure 1 As shown,

[0090] (1) Select limestone aggregates with different particle size ranges and perform CT scanning and 3D modeling reconstruction;

[0091] (101) The aggregate obtained from the quarry is screened and classified into four categories according to the particle size based on the screen aperture size: 2.36-4.75mm, 4.75-9.5mm, 9.5-13.2mm and 13.2-16.0mm.

[0092] (102) At least 25 aggregates of each particle size were selected, washed and fixed in a box. An industrial computer tomography scanner was used to scan the aggregates in the box from front to back to obtain the tomographic projection image of the aggregates.

[0093] In this embodiment, a Phoenix v|tome|xm industrial computed tomography scanner from Germany was used, with a scanning voltage of 200kV, a current of 250μA, a voxel resolution of 40μm for the sample, a magnification of 5x, and a scanning time of 40min.

[0094] (103) The scanned tomographic images were imported into AVIZO software for morphological processing, including filtering, noise reduction, contrast enhancement and other methods. Threshold segmentation was used for binarization to identify the planar regions containing aggregates. The three-dimensional structural models of each aggregate were obtained by using the software's automatic stacking.

[0095] (104) Continue generating a surface mesh covering the 3D aggregate contour in the software. Simplify the aggregate surface mesh by controlling the number of points, surfaces, and edges. The simplification target for a single aggregate surface mesh is 1000 points, 2000 surface faces, and 3000 edges, such as... Figure 2 .

[0096] (2) Import the above aggregate coordinate file into the discrete element software, determine the constitutive model parameters between aggregate and asphalt mortar through uniaxial creep test, and generate virtual samples of asphalt mixture according to the gradation type;

[0097] (201) After simplifying the surface mesh of the aggregate, its contour coordinate information is output to an STL file, and a command is written to import the above STL file into the Discrete Element Method (PFC) software to create an aggregate cluster template. By adjusting the two parameters that affect the simulation accuracy of the aggregate cluster template, namely distance and ratio, a balance between simulation accuracy and simulation efficiency is achieved;

[0098] The distance parameter controls the smoothness of the template, ranging from 0 to 180. Higher values ​​indicate a smoother aggregate surface. The ratio parameter represents the radius ratio of the smallest and largest pebbles within the aggregate template, ranging from 0 to 1. Smaller ratio values ​​allow for more accurate aggregate morphology simulation. After multiple trials, the distance and ratio parameters were determined to be 140 and 0.3, respectively. The number of filling balls for different distance and ratio parameters is compared as follows: Figure 3

[0099] (202) The elastic modulus of the aggregates was obtained by nanoindentation test, and the constitutive model parameters of the contact between the aggregates were determined by selecting a linear stiffness model. The expression is as follows:

[0100] k n =2EL

[0101]

[0102] In the formula, E is the elastic modulus of the aggregate, which is determined to be 55.5 MPa according to the nanoindentation test; L is the distance between the two contact ball elements, L=R1+R2, which is 2 mm in this study; υ is the Poisson's ratio of the aggregate, which is 0.25;

[0103] (203) After constructing the aggregate model, asphalt mortar containing asphalt binder and fine aggregate with a particle size of less than 2.36 mm was prepared in the laboratory and uniaxial creep loading tests were conducted at different temperatures. The creep model parameters were obtained by fitting the test data using Burger's viscoelastic model, and its expression is as follows:

[0104]

[0105] In the formula, E1 is the elastic modulus of the Kelvin portion of the Burger's model; η1 is the viscosity of the Kelvin portion of the Burger's model; E2 is the elastic modulus of the Maxwell portion of the Burger's model; η2 is the viscosity of the Maxwell portion of the Burger's model; and t is the loading time.

[0106] In this example, AC-13 suspended dense gradation mixture was selected, and uniaxial creep loading tests were carried out on asphalt mortar at 40℃, 50℃ and 60℃. The creep parameter fitting is shown in Table 1.

[0107] Table 1 Macroscopic parameters of the Burger's model

[0108]

[0109] Simultaneously, the creep model parameters obtained above are converted into constitutive model parameters for the contact between aggregates and asphalt mortar in discrete element software, as expressed below:

[0110] C mn =η1L

[0111] K mn =E1L

[0112] C kn =η2L

[0113] K kn =E2L

[0114]

[0115]

[0116]

[0117]

[0118] In the formula, C mn K mn C kn and K kn To simulate the parameters of the Burger's model at the normal direction; C ms K ms C ks and K ks υ′ represents the parameters of the Burger's model in the shear direction at the simulation scale; υ′ represents the Poisson's ratio of the asphalt mortar.

[0119] Table 2 summarizes the parameters of the Burger's model after the transformation to the simulated scale.

[0120] Table 2. Burger's model parameters at the simulation scale.

[0121]

[0122] (204) Generate a cuboid box and corresponding six walls based on the size of the rut sample. Calculate the number of four aggregate particles in step (101) according to the required gradation type and put them into the virtual box formed by the rut sample. Scaling and expanding these aggregates in sequence to eliminate the stress caused by particle overlap.

[0123] (205) Generate a specified number of small spheres with a diameter of 1.0 mm to represent asphalt mortar in the remaining space of the rutting sample box, and delete a fixed number of small spheres according to the required porosity of the sample, thus forming a virtual rutting sample of asphalt mixture, such as... Figure 4 .

[0124] (3) Apply load to virtual sample, identify contact between aggregate and asphalt mortar, and classify the contact normal and normal contact force of the above contact according to the angle between the contact vector and the coordinate axis.

[0125] (301) After creating the virtual sample, remove the walls above the sample while retaining the walls on the sides and bottom, and generate a loading plate with dimensions of 50mm × 50mm directly above the sample, such as... Figure 5 ;

[0126] (302) Assign the contact constitutive model parameters to each component, and convert the moving loading time specified in the laboratory rutting test into static loading and determine the equivalent loading time, as shown in the following expression:

[0127]

[0128]

[0129] In the formula, l is the total contact length between the rubber wheel and the sample; m is the mass of the rubber wheel, which is 78 kg; g = 9.8 m / s² 2 , where is the acceleration due to gravity; p is the contact pressure, equal to 0.7 MPa; d = 0.05 m, which is the width of the rubber wheel; t is the cumulative loading time.

[0130] The constitutive model parameters obtained in step (203) were scaled by a factor of 10,000 using the time-temperature equivalence principle, and a program was written to input them into the discrete element PFC software;

[0131] (303) Use the FISH language in PFC software to activate the servo system to apply the load, and record the particle contact between aggregates, between asphalt mortars and between aggregates and asphalt mortars during loading;

[0132] (304) After loading is complete, output particle contact information, including normal contact force f. n and the three-dimensional contact normal vector n i Simultaneously, the three-dimensional contact vector is projected onto the XOZ plane to obtain the two-dimensional projection vector n of the contact normal vector. ip ;

[0133] (305) Calculate the angle θ between each contact vector and the X-axis, ranging from 0° to 360°, and sort these vectors from low to high according to the size of the angle.

[0134] (4) Convert the above contact normal vector into a unit vector, calculate the structure tensor of the contact normal, and convert the contact normal and normal contact force vector into polar coordinates, and draw the wind rose diagram at specified angle intervals.

[0135] (401) Convert the three-dimensional contact normal vector obtained in step (304) into a unit vector, and calculate the configuration tensor of the contact normal, which is expressed as:

[0136]

[0137] In the formula, is the structural tensor; N is the total number of contact normal vectors in the asphalt mixture; is the component of the unit vector of the k-th contact normal; It is the tensor product;

[0138] Taking the contact normal as an example, we define the contact normal probability density function E(n), where the average value of E(n) across all directions in the entire spatial domain is 1, i.e.

[0139] ∫ Ω E(n)dΩ=1

[0140] In the formula, Ω is the solid angle corresponding to the entire surface of a unit sphere.

[0141] Therefore, the discrete contact in the above structural tensor calculation equation can be replaced by certain continuous distributions, expressed as:

[0142]

[0143] E(n) can be expanded using the second-order Fourier method, as follows:

[0144]

[0145] In the formula, D ij yes The deviation portion is represented as:

[0146]

[0147] In the formula, δ ij It's the Kronecker function:

[0148]

[0149] (402) The result is a 3×3 symmetric matrix. The three eigenvalues ​​of the construct tensor can be further extracted ( and The degree of anisotropy is measured by a tensor called the generalized octahedral structure tensor Φ.

[0150]

[0151] (403) The two-dimensional projection vector n obtained in step (304) is... ip The Cartesian coordinates were converted to polar coordinates, and wind rose diagrams were drawn at 10° angular intervals. In the rose diagram of the contact normal, the number of each bin was normalized by the total number of contacts, while in the rose diagram of the normal contact force, the number of each bin was normalized by the average normal contact force of all contacts.

[0152] (5) Construct Fourier series functions and fit wind rose diagrams of contact normals and normal contact forces to evaluate the anisotropy of asphalt mixture skeleton contact using the fitted anisotropy parameters and fabric tensors.

[0153] (501) For the two-dimensional wind rose diagram drawn in step (403), the contact normal density function has a two-dimensional similarity in polar coordinates, namely

[0154]

[0155] In the formula, a ijIt is the coefficient of contact normal anisotropy, equal to:

[0156]

[0157] Distribution function E(θ) and These are statistics of the contact normal and mean normal contact force in polar coordinates. Typically, they can be approximated by a second-order Fourier series:

[0158]

[0159]

[0160] In the formula, a is the second-order anisotropy parameter of the contact normal; θ a The principal direction that touches the anisotropy of the normal; a is the average normal contact force of all contacts in the asphalt mixture; n θ is the second-order anisotropy parameter of the normal contact force; f The principal direction of the anisotropy of the normal contact force;

[0161] (502) In the engineering drawing software OriginPro, the distribution function E(θ) and The two-dimensional wind rose diagram drawn in the fitting (403) step is as follows: Figure 6 and Figure 7 The values ​​of a and a can be obtained through fitting. n Evaluate the two-dimensional anisotropy of the contact normal and normal contact forces;

[0162] (503) The three-dimensional anisotropy of the asphalt mixture skeleton contact was evaluated using the generalized octahedral contact normal tensor Φ obtained in step (402), and the results are shown in Table 3. It can be considered that for AC-13 asphalt mixture, the anisotropy of the contact normal and normal contact force increases with increasing temperature.

[0163] Table 3. Parameters for evaluating two-dimensional and three-dimensional anisotropy

[0164]

[0165] The above embodiments are merely illustrative of the technical concept of the present invention and should not be construed as limiting the scope of protection of the present invention. Any modifications made to the technical solutions based on the technical concept proposed in this invention shall fall within the scope of protection of this invention.

Claims

1. A method for evaluating the anisotropy of skeleton contact of asphalt mixture based on discrete elements, characterized by, Includes the following steps: (1) Select limestone aggregates with different particle size ranges and perform CT scanning and 3D modeling reconstruction; (2) Import the above aggregate coordinate file into the discrete element software, determine the constitutive model parameters between aggregate and asphalt mortar through uniaxial creep test, and generate virtual samples of asphalt mixture according to the gradation type; (3) Apply load to virtual sample, identify contact between aggregate and asphalt mortar, and classify the contact normal and normal contact force of the above contact according to the angle between the contact vector and the coordinate axis. (4) Convert the above contact normal vector into a unit vector, calculate the structure tensor of the contact normal, and convert the contact normal and normal contact force vector into polar coordinates, and draw the wind rose diagram at specified angle intervals. (5) Construct Fourier series functions and fit wind rose diagrams of contact normals and normal contact forces to evaluate the anisotropy of asphalt mixture skeleton contact using the fitted anisotropy parameters and fabric tensors.

2. The method for evaluating asphalt mixture skeleton contact anisotropy according to claim 1, characterized by, In step (1), selecting limestone aggregates of different particle size ranges and performing CT scanning and 3D modeling reconstruction includes: (101) The aggregate obtained from the quarry is screened and classified into four categories according to the particle size based on the screen aperture size: 2.36-4.75mm, 4.75-9.5mm, 9.5-13.2mm and 13.2-16.0mm. (102) At least 25 aggregates of each particle size were selected, washed and fixed in a box. An industrial computer tomography scanner was used to scan the aggregates in the box from front to back to obtain the tomographic projection image of the aggregates. (103) The scanned tomographic images were imported into AVIZO software for morphological processing, including filtering, noise reduction, and contrast enhancement. Threshold segmentation was used for binarization to identify the planar regions containing aggregates. The three-dimensional structural models of each aggregate were obtained by using the software's automatic stacking. (104) Continue to generate a surface mesh covering the three-dimensional aggregate profile in the software, and simplify the surface mesh of the aggregate by controlling the number of points, surfaces and edges.

3. The method of claim 1, wherein the asphalt mixture skeleton contact anisotropy is evaluated by a method comprising: In step (2), the aggregate coordinate file is imported into the discrete element method software. The constitutive model parameters between the aggregate and the asphalt mortar are determined through uniaxial creep tests. Virtual samples of asphalt mixtures are generated according to the gradation type, including: (201) After simplifying the surface mesh of the aggregate, output its contour coordinate information to an STL file, and write a command to import the above STL file into the discrete element PFC software to create an aggregate cluster template; by adjusting the two parameters that affect the simulation accuracy of the aggregate cluster template, namely distance and ratio, a balance between simulation accuracy and simulation efficiency is achieved. (202) The elastic modulus of the aggregates was obtained by nanoindentation test, and the constitutive model parameters of the contact between the aggregates were determined by selecting a linear stiffness model. The expression is as follows: k n = 2EL In the formula, E is the elastic modulus of the aggregate, which is determined to be 55.5 MPa according to the nanoindentation test; L is the distance between the two contact ball elements, L = R1 + R2, which is 2 mm; υ is the Poisson's ratio of the aggregate, which is 0.

25. (203) After constructing the aggregate model, asphalt mortar containing asphalt binder and fine aggregate with a particle size of less than 2.36 mm was prepared in the laboratory and uniaxial creep loading tests were conducted at different temperatures. The creep model parameters were obtained by fitting the test data using Burger's viscoelastic model, and its expression is as follows: In the formula, E1 is the elastic modulus of the Kelvin portion of the Burger's model; η1 is the viscosity of the Kelvin portion of the Burger's model; E2 is the elastic modulus of the Maxwell portion of the Burger's model; η2 is the viscosity of the Maxwell portion of the Burger's model; and t is the loading time. Simultaneously, the creep model parameters obtained above are converted into constitutive model parameters for the contact between aggregates and asphalt mortar in discrete element software, as expressed below: C mn = η1L K mn = E1L C kn = η2L K kn = E2L In the formula, C mn K mn C kn and K kn To simulate the parameters of the Burger's model at the normal direction; C ms K ms C ks and K ks The parameters of the Burger's model at the simulation scale are in the shear direction; υ′ is the Poisson's ratio of the asphalt mortar; (204) Generate a cuboid box and corresponding six walls based on the size of the rut sample. Calculate the number of four aggregate particles in step (101) according to the required gradation type and put them into the virtual box formed by the rut sample. Scaling and expanding these aggregates in sequence to eliminate the stress caused by particle overlap. (205) Generate a specified number of small balls with a diameter of 1.0 mm in the remaining space of the rut sample box to represent asphalt mortar, and delete a fixed number of small balls according to the required porosity of the sample to finally form a virtual rut sample of asphalt mixture.

4. The method of claim 1, wherein the asphalt mixture skeleton contact anisotropy is evaluated by a method comprising: In step (3), a load is applied to the virtual sample to identify the contact between the aggregate and the asphalt mortar. The contact normal and normal contact force of the above contact are classified according to the angle between the contact vector and the coordinate axis, including: (301) After creating the virtual sample, remove the wall above the sample while keeping the walls on the sides and bottom, and generate a loading plate with a size of 50mm×50mm directly above the sample. (302) Assign the contact constitutive model parameters to each component, and convert the moving loading time specified in the laboratory rutting test into static loading and determine the equivalent loading time, as shown in the following expression: In the formula, l is the total contact length between the rubber wheel and the sample; m is the mass of the rubber wheel, which is 78 kg; g = 9.8 m / s² 2 , where is the acceleration due to gravity; p is the contact pressure, equal to 0.7 MPa; d = 0.05 m, which is the width of the rubber wheel; t is the cumulative loading time; The constitutive model parameters obtained in step (203) were scaled by a factor of 10,000 using the time-temperature equivalence principle, and a program was written to input them into the discrete element PFC software; (303) Use the FISH language in PFC software to activate the servo system to apply the load, and record the particle contact between aggregates, between asphalt mortars and between aggregates and asphalt mortars during loading; (304) After loading is complete, output particle contact information, including normal contact force f. n and the three-dimensional contact normal vector n i Simultaneously, the three-dimensional contact vector is projected onto the XOZ plane to obtain the two-dimensional projection vector n of the contact normal vector. ip ; (305) Calculate the angle θ between each contact vector and the X-axis, ranging from 0° to 360°, and sort these vectors from low to high according to the size of the angle.

5. The method of claim 1, wherein the asphalt mixture skeleton contact anisotropy is evaluated by a method comprising: In step (4), the contact normal vector is converted into a unit vector, the configuration tensor of the contact normal is calculated, and the contact normal and normal contact force vectors are converted into polar coordinates. The wind rose diagram is then drawn at specified angular intervals, including: (401) Convert the three-dimensional contact normal vector obtained in step (304) into a unit vector, and calculate the configuration tensor of the contact normal, which is expressed as: In the formula, is the structural tensor; N is the total number of contact normal vectors in the asphalt mixture; is the component of the unit vector of the k-th contact normal; It is the tensor product; Taking the contact normal as an example, we define the contact normal probability density function E(n), where the average value of E(n) across all directions in the entire spatial domain is 1, i.e. In the formula, Ω is the solid angle corresponding to the entire surface of a unit sphere; Therefore, the discrete contacts in the above structural tensor calculation equation are replaced by certain continuous distributions, expressed as: E(n) can be expanded using the second-order Fourier method, as follows: where D ij is the bias part, expressed as: where δ ij is the Kronecker function: (402) The result is a 3x3 symmetric matrix; further extraction of the three eigenvalues of the fabric tensor and to measure the degree of anisotropy, called the generalized octahedral fabric tensor Φ: (403) The two-dimensional projection vector n obtained in step (304) is... ip The Cartesian coordinates were converted to polar coordinates, and the wind rose diagrams were drawn at 10° angular intervals. In the rose diagram of the contact normal, the number of each bin was normalized by the total number of contacts, while in the rose diagram of the normal contact force, the number of each bin was normalized by the average normal contact force of all contacts.

6. The method of claim 1, wherein the asphalt mixture skeleton contact anisotropy is evaluated by a method comprising: In step (5), a Fourier series function is constructed and a wind rose diagram of the contact normal and normal contact force is fitted. The anisotropy parameters and fabric tensor obtained from the fitting are used to evaluate the anisotropy of the asphalt mixture skeleton contact, including: (501) For the two-dimensional wind rose diagram drawn in step (403), the contact normal density function has a two-dimensional similarity in polar coordinates, namely where a ij is the coefficient of contact normal anisotropy, equal to: The distribution functions E(θ) and are statistical quantities of the contact normal and the mean normal contact force in polar coordinates; they are approximated by a second order Fourier series: In the formula, a is the second-order anisotropy parameter of the contact normal; θ a The principal direction that touches the anisotropy of the normal; a is the average normal contact force of all contacts in the asphalt mixture; n θ is the second-order anisotropy parameter of the normal contact force; f The principal direction of the anisotropy of the normal contact force; (502) In the engineering drawing software Origin Pro, the distribution function E(θ) and The two-dimensional wind rose diagram drawn in step (403) is fitted, and a and a are obtained through fitting. n Evaluate the two-dimensional anisotropy of the contact normal and normal contact forces; (503) The three-dimensional anisotropy of the asphalt mixture skeleton contact is evaluated using the generalized octahedral contact normal tensor Φ obtained in step (402).

Images