Figure 16 shows the results of fitting the BPN value and MPD of the two mixtures. The results show that the BPN value and MPD also have a good linear relationship. The goodness of fit of the two mixtures is above 0.84. The average section depth of MPD decreases with the decrease of the BPN value. The analysis shows that as the number of aeolian sand increases, the texture depth of the specimen surface gradually decreases, and the surface texture changes passively. At this time, the MPD decreases accordingly.

**Figure 16.** Relationship between BPN and MPD.

#### **4. Conclusions**

This paper studies the influence of aeolian sand on the anti-sliding performance and texture index of asphalt mixture, which is a study based on the macro perspective of the anti-sliding performance changes caused by different amounts of sand and asphalt mixture types.

**Author Contributions:** Conceptualization, J.P. and H.Z.; methodology, J.P.; literature survey, J.P. and H.Z.; validation, J.P., H.Z. and G.L.; formal analysis, J.P. and H.Z.; investigation, J.P., H.Z. and Y.W.; resources, Y.W. and G.L.; data curation, J.P. and H.Z.; writing—original draft preparation, J.P.; writing—review and editing, G.L. and H.Z.; supervision, Y.W. and G.L.; project administration, Y.W. and G.L.; funding acquisition, Y.W. and G.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Fund of China, grant numbers 51768062 and 52068061.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No new data were created or analyzed in this study, data sharing is not applicable to this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Anti-Skid Characteristics of Asphalt Pavement Based on Partial Tire Aquaplane Conditions**

**Miao Yu 1,2,\*, Yao Kong 3, Zhanping You 4,\*, Jue Li 2, Liming Yang 5 and Lingyun Kong 2**

**Abstract:** This study presented a finite element model of radial tire–asphalt pavement interaction using ABAQUS 6.14 software to investigate the skid resistance properties of asphalt pavement under partial tire aquaplane conditions. Firstly, the pavement profile datum acquired by laser scanning were imported to Finite Element Analysis (FEA) software to conduct the pavement modeling. Secondly, a steady state rolling analysis of a tire on three types of asphalt pavements under drying conditions was carried out. Variation laws of the friction coefficient of the radial tire on different pavements with different pavement textures, tire pressures, and loads on the tire were examined. Subsequently, calculation results of the steady state rolling analysis were transmitted to dynamic explicit analysis, and an aquaplane model of a radial tire on asphalt pavements was built by inputting the flow Euler grids. The tire–pavement adhesive characteristics under partial aquaplane conditions are discussed regarding the aquaplane model. Influences of the thickness of water film, the texture of asphalt pavement, and the rolling speed of the tire on the vertical pavement-tire contact force are analyzed. It is found that the vertical contact force between open graded friction course (OGFC) pavement and tire is the highest, followed by stone mastic asphalt (SMA) pavement and dense graded asphalt concrete (AC) pavement surface. The vertical contact force between tire and pavement will be greatly reduced, even with increasing speed or water film thickness. As tire speed increases from 70 km/h to 130 km/h, the tire–pavement contact force is reduced by about 25%. Moreover, when the thickness of water film increases from 0 (dry condition) to 4 mm and then to 12 mm, the vertical contact force reduced 50% and 15%, respectively, compared with under the dry contact condition. This study provided a key theoretical reference for safe driving on wet pavements.

**Keywords:** skid resistance; tire-pavement friction; pavement textures; partial tire aquaplane conditions; water film; vertical contact force

#### **1. Introduction**

The anti-skid performance of pavement is one of the major influencing factors of traffic accidents. Inadequate skid-resistance may decrease the braking efficiency of vehicles, contributing to accidents, including collisions [1,2]. Whether exhibiting the pure skidding behavior of lock slip or the usual simultaneous rolling-skidding state of vehicles, the tire– pavement friction behavior is the control factor of the braking effect since the braking force of the tire on the pavement is significantly higher than the rolling resistance in the braking process [3,4]. Under severe weather conditions like heavy rain, the pavement surface is covered by a layer of water film, as the rainwater cannot be discharged quickly

**Citation:** Yu, M.; Kong, Y.; You, Z.; Li, J.; Yang, L.; Kong, L. Anti-Skid Characteristics of Asphalt Pavement Based on Partial Tire Aquaplane Conditions. *Materials* **2022**, *15*, 4976. https://doi.org/10.3390/ma15144976

Academic Editor: Francesco Canestrari

Received: 19 June 2022 Accepted: 15 July 2022 Published: 17 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

enough. When vehicles drive on the pavement covered with water film, the tires are lifted by water flow and cannot contact the pavement sufficiently. The adhesive characteristics and braking characteristics of vehicles decrease significantly compared with those on dry pavement, thus greatly increasing the risk of traffic accidents [5]. In addition, traffic accidents induced by weather conditions not only have a high frequency of occurrence, but also are accompanied by secondary accidents [6,7]. Weather-related traffic accidents have attracted much attention from the traffic safety industry [8,9]. Once the tire aquaplane phenomenon occurs, adhesive characteristics of tires and the anti-skid characteristics of pavement under wet conditions is crucial to ensuring driving safety [10,11]. However, anti-skid performance of asphalt pavement under tire aquaplane conditions is mainly discussed through simulation analysis due to the great differences between test simulation and practical vehicle-pavement coupling friction conditions.

Ong et al. [12] built a finite element model to evaluate anti-skid performance of trucks on wet pavement. The model covers basic structural mechanics and fluid mechanic theory, as well as tire contact and tire–fluid interaction. The model is used to analyze the effects of speed, tire load, road thickness, and water film thickness on the skid resistance of trucks under zero load and full load. With consideration to the influences of the geometric design of pavement, tire surface design, and various working conditions of tires, a finite element model of tire–water–pavement interaction was established by Tang Tianci et al. [13] to evaluate the slipping performance of vehicles under different rainfall intensities. The real pavement surface texture and the internal microstructure of the asphalt mixture were integrated into the pavement grid by CT scanning. The established model could not only simulate water flow on a pavement surface, but also capture the process of water penetration into the porous structure of pavement, thus enabling the simulation of the tire–water–pavement interaction effectively. Furthermore, Anupam [14] also conducted an investigation on the tire–pavement heat transfer by modeling the coupling contact between tire and asphalt pavement. The textured pavement surface submodel were modeled by the combination of CT scanning and the commercial software Simpleware, which was incorporated for the post-process of surface images. Feng et al. [15] constructed finite element models of tires with smooth, longitudinal grooves and tires with complete texture in ABAQUS by using the three-dimensional solid modeling software CATIA. These models were used for water-skidding simulations of tires under static inflation loading, drive braking, free rolling, and transient water-skidding states, respectively, which further verified the validity of the CEL method in studying water-skidding behaviors of tires. Zhu et al. [16] built an aquaplane model of rib tires on asphalt pavement by using ABAQUS. The solid model of asphalt pavement was established by CT scanning, accompanied by the coupled tire model and fluid model. Furthermore, the effect of tire pressure, thickness of water film, and macrotexture of asphalt pavement on a tire aquaplane were analyzed. Zhu et al. [17] built an aquaplane model of rib tires on runways in airports by using the CEL algorithm with consideration to the macrotexture of asphalt pavement, and they used the model to discuss influences of slip ratio, flow distribution, and tire pressure on a tire aquaplane. Numerical simulation results were verified by calculating the braking distance and trailer test.

To sum up, there are various influencing factors of tire–water–pavement coupling friction. For simplification of calculation in road engineering, pavement is usually simplified as a rigid plane [18]. That is, the investigation of wet friction performance of pavement is insufficient in comparison. With further deepening studies, some scholars have attempted to study the tire–pavement coupling friction behaviors in humid environments through CT scanning and have achieved good evaluation effects [19,20]. However, these studies often have high test costs and require much time for calculation, which brings some objective limitations against further investigation.

Based on the above brief review of the research status of tire–pavement contact behaviors on tire aquaplane conditions, it can be inferred that the research on pavement skid resistance, especially when tires rolling on partial aquaplane conditions is insufficient. In

addition, the texture characteristics are acquired mainly by CT scanning, which also incurs high cost. Therefore, the inadequate investigation of tire rolling on wet pavement surfaces with different thicknesses of water film has prevented further study on pavement frictional mechanism under the wet condition of pavement design, which also plays a key role in traffic safety.

The objective of this study is to investigate the frictional properties of asphalt pavement on partial tire aquaplane conditions at a relatively low cost. The pavement texture characteristics can be acquired by laser scanning, and then the textured pavement can also be modeled in ABAQUS by data processing. Variation laws of the friction coefficient of the radial tire on various pavements with the associated pavement textures, tire pressures, and loads on the tire are further investigated. Moreover, the tire–pavement contact characteristics under partial aquaplane conditions are discussed in explicit analysis mode. Influences of the thickness of water film, the texture of asphalt pavement, and the rolling speed of the tire on the vertical pavement-tire contact force are analyzed. It seeks a relatively good computer analysis target at a low computational cost and provides application references as well as theoretical references to the anti-skid performance of the wet pavement.

#### **2. Methodology**

Based on the finite element model of radial tire–asphalt pavement interaction, which was built with ABAQUS, a steady state rolling analysis of a tire on three drying asphalt pavements was carried out first to explore variation laws of the friction coefficient of radial tires with pavement texture, tire pressure, and loads on tire. Second, calculation results of the steady state rolling analysis were transmitted to a dynamic explicit analysis, and an aquaplane model of a radial tire on asphalt pavements was built by inputting the flow Euler grids. The model is used to analyze the adhesion characteristics of tire and pavement under partial aquaplane conditions and explore the effects of water film thickness, asphalt pavement texture, and tire movement speed on the vertical contact force between pavement and tire, so as to provide some theoretical guidance for the safety of vehicle driving on rainy days.

#### *2.1. Construction of the Tire–Pavement Contact Model*

#### 2.1.1. Tire Modeling and Assembly with Pavement

In this study, the tire–pavement contact simulation analysis was carried out using ABAQUS finite element simulation software. The tire model used the 175SR14 radial tire. Initially, the axisymmetric model is constructed according to the geometry and material properties of a radial tire. Subsequently, the partial 3D model and an integrated threedimensional model of the tire can be generated in turn in ABAQUS, as shown in Figure 1.

**Figure 1.** Tire modeling. (**a**) Axisymmetric model; (**b**) partial 3D model; (**c**) full 3D model.

Based on the Input file of the completed 3D tire model, the keyword \*INCLUDE, INPUT was adopted for importing the Input file of the asphalt pavement model. Moreover, the tire–pavement contact behaviors were defined, accompanied with the set of loading and boundary conditions, to realize the tire–pavement assembly, as illustrated in Figure 2. More details on modeling and assembly can be found in Reference [21]. In addition, the construction steps of the pavement texture model are illustrated as follows.

**Figure 2.** Tire–pavement assembly.

#### 2.1.2. Acquisition of Pavement Texture Properties

The laser profile scanner was assembled to acquire relative altitude information on surfaces of asphalt concrete specimens. This laser scanner was assembled with the crossskidding table and double-shaft connected control cabinet, as shown in Figure 3a. To acquire sufficient texture properties of asphalt pavement, a 100 mm × 100 mm rectangular zone of each specimen was scanned, with a 1 mm points-acquisition interval in the scanning area, as shown in Figure 3b,c. That is, 10,000 effective altitude data were collected in each scanning region.

**Figure 3.** *Cont*.

**Figure 3.** Laser scanning of asphalt pavement profile. (**a**) Basic configuration of laser scanner; (**b**) scanning path; (**c**) scanning region of texture.

The laser displacement sensor moves continuously along the path in Figure 3b, according to the programming of the control cabinet. The output data of the supporting software are continuous two-dimensional altitude data, and the continuous altitude data were, thereafter, divided into different contour lines. Since the standard scanning range of the laser scanning head is ±10 mm, the data were uniformly displayed as 12.5 mm. Therefore, the singular points should be deleted and interpolated prior to data processing. Following that, in order to differentiate the contour lines, the transverse length of the scanning path was extended to get the significant height difference of the altitude, and the boundary of the contour line could be easily verified according to the data mutation. Moreover, altitude data of the even-order contour lines were reverted around the central point and then arranged according to the relative positional relations of the contour lines. Thus, 100 contour altitude lines with an interval of 1 mm along the same direction were obtained.

To eliminate the uneven phenomena of the asphalt pavement surface originating from specimen compaction, linear fitting was conducted to each contour line. Moreover, 100 altitude data of each contour line were arranged in order from small to high, and the absolute value of the median was set as the elevation zero point. Subsequently, the whole contour line was moved vertically to eliminate influences of average altitude differences among specimens. Following that, the processed texture data were input into the MAT-LAB software to reconstruct three-dimensional profile maps of the three types of asphalt pavement, namely, Asphalt Concrete (AC), Stone Matrix Asphalt (SMA), and Open Graded Friction Course (OGFC), as shown in Figure 4.

#### 2.1.3. Modeling of Asphalt Pavement with Profile Textures

The collected three-dimensional coordinates of asphalt pavement were transformed into a three-dimensional grid model through Python programming to generate the INP file of pavement for further simulation analysis in ABAQUS.

In general, the reconstructed three-dimensional geometric models of asphalt pavement can be categorized as wireframe models, curved shell models, and solid models. Only the undulating state of pavement surface is needed in the modeling of curved shell; that is, the whole rigid surface can be constituted by a series of triangular or quadrilateral rigid slices so as to realize the accurate characterization of road surface texture. Hence, this type of shell model is adopted in this research to build the 3D model of asphalt pavement with the following steps.

**Figure 4.** Three—dimensional profile maps of different types of asphalt pavement (Unit: mm). (**a**) AC pavement; (**b**) SMA pavement; (**c**) OGFC pavement.

(1) Pavement filtering and smoothing

Considering that a large number of sharp angles can be generated from the curved shell modeling based on the data of pavement texture, significant stress concentration often occurs in the areas that contact those sharp peaks of the pavement surface. This phenomenon always leads to the penetration of the road model to the tire model exceeding the tolerance range, thus triggering the termination of model calculation. In addition, the tire–pavement contact still cannot converge if the meshing of the pavement model is too fine. As a result, it was necessary to implement pavement filtering and smoothing processing first to eliminate sharp points on the pavement surface. Next, the pavement was amplified to some extent along the plane *xoy* to make the grid size of the pavement larger than that of the tire. Furthermore, mean filtering and Grid data function were applied in combination to conduct the smoothing process to the texture data, as illustrated in Figure 5.

(2) Extension and amplification of pavement

Compared with the dimensions of the tire model, the scanning region of the asphalt pavement was too small to develop normal contact with the tire. Due to the self-similarity and self-affinity of asphalt pavement, the pavement model could be expanded both in longitudinal and lateral directions. To ensure a smooth transition of altitude between the original pavement and the expanded pavement, the altitude data of the original pavement were expanded asymmetrically along the pavement boundary. Subsequently, the original pavement was expanded by three times along the longitudinal and transverse directions, respectively. Meanwhile, local coordinate systems were defined through the keyword of \*SYSTEM, and the center of the expanded pavement was moved to the place below the tire.

**Figure 5.** Grid data diagram of pavement profile in 150 × 150 area.

The dimensions of the pavement grid model after extension were 300 mm × 300 mm, which covered 150 × 150 altitude data points, with 2 mm precision. In contrast, the size of the texture unit of the tire was about 2.5 mm. Since the grid size of the pavement was smaller than the grid size of the tire, the tire–pavement contact analysis could not reach convergence. Hence, it was necessary to expand the pavement coordinates by 1.5 times along *x*, *y*, and *z* directions under the premise of an unchanged number of altitude data points. Therefore, the size of the pavement grid was set to 3 mm, and the whole size of the pavement model became 450 mm × 450 mm, as shown in Figure 6.

(3) Generation of the grid model

Altitude data in Microsoft EXCEL were extracted by Python language programming and then converted into 3D coordinates. Four adjacent nodes formed a quadrilateral element, and corresponding relations between node number and element number were established, as presented in Figure 7. All nodes and elements were written into the keyword \*NODE of nodes in turn, and the keyword \*ELEMENT for element was defined. The element type was set as R3D4. Finally, the INP format file of the finite element model of pavement was generated, which contained 22,500 nodes and 22,200 elements. The finite element grid models of a 100 mm × 100 mm pavement after mean filtering and smoothing before extension and amplification are shown in Figure 8.

**Figure 6.** Pavement grid setting after the extension and the amplification processing.

**Figure 7.** Schematic diagram of node number and element number.

**Figure 8.** Finite element grid models of asphalt pavements. (**a**) AC pavement; (**b**) SMA pavement; (**c**) OGFC pavement.

Since the rolling speed of the tire cannot be defined directly in steady-state rolling analysis, it is necessary to determine the rolling angular velocity and linear velocity of the tire through the keywords \*TRANSPORT VELOCITY and \*MOTION. Specifically, \*TRANSPORT VELOCITY is used to define the angular velocity of tire material running through the model grid. \*MOTION is used to define the kinematic velocity of the moving reference coordinates, which is the translational velocity of the tire. As the tire rolls freely, the linear velocity of the tire (V) is just equal to the product of the rolling angular velocity (ω0) and rolling radius (r). In addition, when the rolling angular velocity of the tire (ω) is smaller than ω0, the tire is in the braking state. By contrast, while ω > ω0, the tire is in the driving state. Hence, the rolling contact analyses of the tire under braking, free rolling, and driving states can be realized by defining ω and V. It should be noted that the rolling radius of the tire is not the radius of the tire after inflation. In other words, static trace analysis is needed before the rolling analysis. Following that, the rolling radius of the tire can be calculated according to the inflating radius of the tire and the radial deformation under the corresponding load. More details about the establishment of the tire–pavement frictional contact simulation model, such as the tire modeling and the definition of tire–pavement contact, can be found at the reference [19].

2.2.2. Influencing Factors of Anti-Skid Performance of Asphalt Pavements

(1) Effects of pavement texture on the frictional force of pavement

The tire pressure was set to 2.5 bar, and a load of 3900 N was applied onto the tire. On this basis, the rolling radii of the tire on AC, SMA, and OGFC were calculated to be 297.91 mm, 297.77 mm, and 297.21 mm, respectively. The steady state rolling analysis when the velocity and slip rate of the tire was 30 km/h and 20% was carried out by defining the rotating angular velocity of the tire and the pavement-to-tire average velocity. The keywords \*TRANSPORT VELOCITY valued 22.4 rad/s and \*MOTION,TYPE = VELOCITY valued 8.3334 m/s.

Figure 9 shows the distribution of the shear force on the rolling direction of the tire on different types of pavements, namely, AC, SMA, and OGFC. The shear force on the tire– pavement interface shows significant differences among the different types of pavement. Specifically, the maximum frictional force that the OGFC provides to the tire element is higher than those of SMA and AC. The frictional shear force of AC pavement is the smallest, at nearly 10 times lower than that of OGFC. Additionally, regions where the tire generates frictional force are closely related to texture features of asphalt pavements. When the pavement texture is relatively flat, frictional force shows a relatively uniform distribution, and it is generated in most regions at the tire surface and tire shoulder in the contact area. Moreover, frictional force at the area with protruding pavement texture is relatively high, indicating that frictional force when the tire rolls on asphalt pavements is mainly provided by protruding textures of the asphalt surface. Furthermore, as pavement texture becomes relatively rough, frictional forces usually show the trend of the dispersing distribution, as illustrated in Figure 9b,c. Since a few coarse aggregates in OGFC are too protruding, frictional forces of the tire at these nodes are negative, thus resulting in the direction of friction force being opposite to the direction of resultant friction force in these areas.

(2) Effects of tire pressure on friction coefficient

In the process of vehicle braking, controlling the slip rate of the tire is conducive to increasing the braking efficiency. The slip rate of the vehicle tire is generally controlled in the range between 15% and 20% to avoid locking of tires during emergency braking and thereby increase braking efficiency. In this study, the variation laws of friction coefficient with tire pressure and loads under free rolling and braking conditions were analyzed by defining the rotating angular velocity of the tire and translational velocity of the pavement. The testing conditions in this simulation were set as follows: tire pressure was set in four ranges, that is, 1.5 bar, 2.0 bar, 2.5 bar, and 3.0 bar, and the vertical loads and translational velocity of pavement were 3900 N and 30 km/h, accompanied by a constant 20% slip ratio. Furthermore, the shear force and vertical reaction force on the tire were output through the keyword \*CONTACT PRINT from the DAT file, thus obtaining the friction coefficient of the tire under different motion states.

It can be seen from Figure 10 that under the free rolling state, all the friction coefficients of the tire on the three types of pavement decrease to some extent with the increase of tire pressure. Additionally, the SMA and OGFC show a significantly greater decreasing trend than the AC, which can be attributed to the rapid reduction of the actual contact area between the tire and the pavement with a rough profile. By contrast, the effects of tire pressure on the friction coefficient under the braking state at a slip rate of 20% are opposite those under the free rolling state. In other words, the friction coefficients of the tire on all three pavements increase with the increase of tire pressure. Moreover, it is found that the friction coefficients of the OGFC and SMA are significantly higher than that of AC pavement under braking conditions. This also sheds light on the fact that pavement texture influences anti-skid performance of pavement significantly.

**Figure 9.** Cloud chart of shear force distribution on different pavements. (**a**) AC Pavement (**b**) SMA Pavement (**c**) OGFC Pavement.

**Figure 10.** Variation laws of friction coefficient with tire pressure. (**a**) Tire free rolling state; (**b**) tire braking state.

#### (3) Effects of loads on the friction coefficient

In this analyzing phase, the simulation conditions were set as follows. Tire pressure and the translational speed of the road relative to the tire were set at 2.5 bar and 30 km/h. Different loads of 1700 N, 2200 N, 2700 N, and 3200 N were applied onto the tire under the free rolling state and braking state with a slip rate of 20%. The dynamic friction coefficients of the tire on different pavements were obtained.

In Figure 11, the variation trend of friction coefficients of the tire on three pavements, both under the free rolling state and under the braking state with a slip rate of 20%, shows positive trends related to the rising loads. In other words, given a certain range of loads, the tire–pavement contact is better with the increase of loads, and, furthermore, the rising trend of the frictional force of the tire provided by pavement is higher than that of the vertical reaction force. In addition, it can be seen from a comparison of the friction coefficients of different pavements that OGFC and SMA have similar friction coefficients, while AC shows a relatively lower friction coefficient, especially under the braking state. Such a difference is mainly due to the different degrees of surface roughness and the corresponding change of actual contact area.

**Figure 11.** Variation laws of friction coefficients with loads. (**a**) Tire free rolling state; (**b**) tire braking state.

#### *2.3. Finite Element Modeling of the Tire Aquaplane*

#### 2.3.1. Principle of CEL Method

In ABAQUS, the fluid-solid coupling issue is mainly solved through the fluid-solid coupling algorithm (CEL), smoothed particle hydrodynamics (SPH), and computational fluid dynamics (CFD). The Euler-Language coupling (CEL) method is adopted to solve the tire aquaplane problem. In CEL, the tire and pavement are defined by Lagrange elements, while the water body is defined by Euler elements, as presented in Figure 12. CEL method combines advantages of the Lagrange method and Euler method. It not only solves the great deformation of fluid grids, but also identifies coupling boundaries effectively. During modeling, the tire grids and fluid grids overlapped to some extent. On fluid free boundaries, pressure on the fluid surface was applied onto the tire surface as intensity of pressure. Meanwhile, the volume fraction and velocity boundary of the fluid grids were updated according to the displacement and velocity of the tire on the coupling surface.

#### 2.3.2. Euler Grid Modeling Based on Flow Model

#### (1) Selection of tire rolling model and flow model

The tire aquaplane model can be categorized as "tire rolling model" and "flow model" according to different selections of reference coordinate systems. It can be seen from Figure 13 that the tire rolling model constrains all degrees of freedom (DOFs) of the pavement model, and it applies the rolling velocity and translational velocity onto the

tire to make it roll on the water film model. The flow model constrains the DOF of the tire in the translational direction, and it applies rotating angular velocity onto the tire only. Meanwhile, translational velocities opposite the advancing direction of the tire were applied onto the pavement model and flow model simultaneously to simulate a tire driving through the water film on pavements. The "tire rolling model" and "flow model" have been compared by some scholars [16]. It is found that the tire rolling model has advantages in simulating flow traces, and it can simulate water traces completely. However, the quantity scale of grids can be decreased significantly during modeling and thereby increases the solving efficiency. Since the tire rolling model and flow model show almost equal calculation accuracies and the flow trace analysis is not the research focus of this paper, the flow model was adopted for the modeling analysis of the tire aquaplane.

**Figure 12.** Tire aquaplane model on asphalt pavements.

**Figure 13.** Schematic diagram of different modeling techniques. (**a**) Tire rolling model; (**b**) flow model.

#### (2) Flow state equation

In this study, the relations of fluid pressure with density and internal energy of fluid were described by the Mie–Gruneisen state equation. It is assumed that the fluid pressure (p) is linear with the internal energy per unit mass (Em):

$$\mathbf{p} = \mathbf{f}(\boldsymbol{\rho}) + \mathbf{g}(\boldsymbol{\rho})\mathbf{E}\_{\mathrm{m}} \tag{1}$$

where f(ρ) and g(ρ) are determined by the selected state equation. In the Mie–Gruneisen state equation, values of f(ρ) and g(ρ) are determined as follows:

$$\mathbf{f}(\boldsymbol{\rho}) = \frac{\mathbf{p}\_0 \mathbf{c}\_0^2 \boldsymbol{\eta}}{\left(1 - \mathbf{s}\boldsymbol{\eta}\right)^2} (1 - \frac{\Gamma\_0 \boldsymbol{\eta}}{2}) \tag{2}$$

$$\mathbf{g}(\rho) = \Gamma\_0 \rho\_0 \tag{3}$$

where *η* = 1 − *ρ*0/*ρ*, *ρ*0 is the initial density of the water body, and *ρ* is the density of water bodies after impacts. c0 is the propagation velocity of acoustic waves in water, and s is a material constant. c0 and s determine the linear relationship between the impact speed of water flow (Us) and particle velocity of water flow (Up): U*s* = c0 + sUp. Г0 is the material constant. The flow parameters in this study are listed in Table 1.

**Table 1.** Parameters of flow materials.

#### (3) Construction of flow grid models

The Euler grid model established in this study was divided into air zone and flow zone. The overall dimensions were 400 mm in *x* direction, 400 mm in *y* direction, and 100 mm in *z* direction, as shown in Figure 14a. To set boundary conditions of the Euler grid model, the flow zone was divided into side, bottom, and inlet of water pump, as shown in Figure 14b,c. The Inlet dimension was 4 mm in *x* direction and 400 mm in *y* direction. Size in the *z* direction the dimension was determined according to the preset thickness of the water film in the model. Furthermore, regions overlapping the tire grid that might be flowed through by water were refined, and the minimum side of the minimum grid was 4 mm. The Euler elements EC3D8R were used as the flow model. The whole flow model contained 110,594 nodes and 100,799 elements.

**Figure 14.** Fluid finite element model. (**a**) Euler grid model; (**b**) side and bottom of flow; (**c**) inlet of water pump (inlet).

2.3.3. Tire Aquaplane Analysis Based on ABAQUS/Standard and Explicit

To reduce the calculation expense in tire aquaplane analysis, an implicit expression was carried out first by using the Standard solver. Next, a dynamic explicit analysis was conducted. Steady state calculation results were transmitted to the explicit analysis by the keyword \*IMPORT, and finally, kinetic analysis in the tire aquaplane process was performed through the Explicit solver.

(1) Implicit analysis

First, a complete 3D radial tire was built in the implicit analysis. The pavement grid files were input into the tire model through the keyword \*INCLUDE, and the tire–pavement contact properties were defined. Static analysis and steady state rolling analysis of the tire on asphalt pavements were carried out. During steady state rolling analysis of the tire, only the meshing of the tire–ground contact area was refined to increase computational efficiency. However, the uniform meshing of the tire is needed in the explicit analysis. Hence, one grid was generated every 4◦ during construction of a complete 3D radial tire model using the keyword \*SYMMETRIC MODEL GENERATION, and a total of 90 tire grids were generated. Since the tire did not achieve real rolling and only materials flow in grids during steady state rolling analysis, a 450 mm × 450 mm asphalt pavement model was constructed. However, the real relative movements were made by the tire and pavement in the explicit analysis. Considering that the asphalt pavements previously built for Standard analysis were too short in size, they had to be extended to 1800 mm × 450 mm, as shown in Figure 15.

**Figure 15.** Modeling process of tire aquaplane.

(2) Explicit analysis

In the explicit analysis, calculation results and pavement grid files of steady state rolling analysis were input from Standard analysis. The concentrated mass and rotational inertia as well as initial velocity of the tire were defined. Meanwhile, the boundary conditions and amplitude of loads were defined by the keyword \*AMPLITUDE. Later, Euler grid files were input through the keyword \*INCLUDE, and material properties of fluid were defined. Additionally, a vertically downward gravitational field was applied to the whole Euler model. When defining the initial state of the fluid, the volume fraction of element set INLET of the water pump was set at 1, while the volume fraction of other regions was set at 0. In other words, only grids of INLET were filled with water completely, whereas other Euler grids were filled with air, as presented in Figure 16, where the red zone is fluid and the blue zone is air.

**Figure 16.** Initial state of tire aquaplane.

In the explicit analysis, the contact relation is defined by the universal contact algorithm. The tire–water contact and tire–pavement contact were built by the keyword \*CONTACT INCLUSIONS. Later, tire–water contact constraint and water–pavement contact constraint were eliminated using the keyword \*CONTACT EXCLUSIONS. In this way, the universal contact properties between the tire and pavement were defined. Subsequently, the contact force of the tire surface which was provided by the tire was output in the historical output through the keyword \*CONTACT OUTPUT.

2.3.4. Validity Verification of Tire Aquaplane Finite Element Model

#### (1) Flow trace verification

In this simulation phase, first, the tire pressure and thickness of water film were set at 2.5 bar and 12 mm, and a vertical load of 3900 N was applied onto the tire. Next, the tire rolled on OGFC was defined at a constant velocity of 70 km/h. The time of analytical step was determined to be 0.03 s. As the water flow impacted the tire, the flow distribution pattern in the Euler grids was disclosed through the field output variable EVF, as presented in Figure 17.

**Figure 17.** Flow trace in the tire aquaplane process. (**a**) t = 0.0060 s; (**b**) t = 0.0105 s; (**c**) t = 0.0165 s; (**d**) t = 0.0225 s.

Figure 17a presents the flow trace distribution as the time history was 0.0060 s. At this time, the tire began to contact the water flow, and the water first flew into longitudinal grooves on the tire surface. The analysis time (t) in Figure 17b was 0.0105 s, when there was no complete contact between the tire and pavement due to the support force of water flow. Specifically, the tire shoulder was in relatively full contact with the pavement, whereas there were some flows between the tire tread and pavement. In other words, loads on the tire were mainly undertaken by the tire shoulder. As the time duration reached 0.0165 s, it can be seen from Figure 17c that the tire surface began to contact the pavement gradually. The longitudinal grooves on the tire surface could not discharge the water between the tire and pavement anymore. Furthermore, the water flow in the longitudinal texture of the tire was dispersed due to the uneven pavement surface. The time duration in Figure 17d is 0.0225 s. In this phase, water flows occupied the whole Euler grid bottom and furthermore, began to splash behind the tire by inertia after tire rolling. Water flows generated the maximum support force upon the squeezing effect, which was kept unchanged. Hence, the aquaplane model in this study could display flow trace distributions in each stage of the tire aquaplane process, which was consistent with practical tire aquaplane conditions.

(2) Vertical contact force

Two tire velocities were set, while other boundary conditions were fixed. The variations of the vertical contact force between the tire and pavement under different tire speeds with time in the aquaplane process were output by defining the historical output, with the fixed boundary conditions, as shown in Figure 18. Before contact with water flows, the tire–pavement contact force fluctuated around 3900 N, which was determined by the damping structure of the tire. With the continuous contact between the tire and water, the lifting effect of the water on the tire led to the gradual decrease of the tire–road contact force. After the tire entered into the wet asphalt pavements, it reached the stress balance state. As the velocity reached 70 km/h, the tire–pavement vertical contact force reached the balance at about 0.02 s, and, finally, it was stabilized at about 1000 N. When the velocity was 150 km/h, the vertical contact force reached the balance at about 0.015 s, and, finally, it approached 0, indicating that the tire was lifted up completely by water flows. In other words, a complete aquaplane phenomenon occurred.

**Figure 18.** Vertical contact force in the tire aquaplane process.

#### **3. Results and Analysis**

Existing studies on tire aquaplane mainly concentrate on complete aquaplane. Experimental studies and numerical models have been conducted by scholars to solve the critical aquaplane velocity and the affecting factors under different conditions [20–22]. Nevertheless, drivers are often more cautious when driving on wet pavements. In other words, the vehicle is usually driven at a low speed, and the tires are in a partial aquaplane

state [23,24], which also may often trigger low adhesive force and traffic accidents. Since the tire–pavement friction behaviors under the partial aquaplane state are closely related with vertical contact force, the effects of pavement types, tire speed, and water film thickness on the vertical contact force of the pavement are investigated as follows.

#### *3.1. Effects of Asphalt Pavement Types on Partial Aquaplane Performanc*

With the same testing conditions applied in the above analysis, the variation trends of vertical contact force from three types of pavement, namely, AC, SMA, and OGFC, are presented in Figure 19, based on history output. All the simulations in this phase were conducted under the same tire speed, 70 km/h.

**Figure 19.** Effects of pavement type on vertical contact force.

According to the variation trend of vertical contact force, there was little difference of vertical contact force among the pavements before the tire–flow contact. However, vertical contact force decreased sharply as the tire entered the wet pavements. Specifically, vertical contact force on AC declined the most, followed by that of SMA and OGFC, which were relatively stable. As the tire reached the equilibrium state of force, the difference of vertical contact force became more significant. OGFC provides the tire with distinct vertical contact force, which is superior to SMA and AC. This difference can be attributed to water discharge capacity, which benefits from the rough surface supplied by sufficient textures, especially on OGFC and SMA pavement.

#### *3.2. Effects of Tire Rolling Velocity on Partial Aquaplane Performance*

The tire pressure, loads, and thickness of water film were set at 2.5 bar, 3900 N, and 12 mm, respectively. With all other parameters being constant, the tire speed was set at 70 km/h, 90 km/h, 110 km/h, and 130 km/h, for the investigation of the effect of tire rolling velocity on partial aquaplane properties. Based on historical output, the variation trends of tire–pavement vertical contact force on AC, SMA, and OGFC in the tire partial aquaplane process with velocity are presented in Figure 20.

According to variation trends of vertical contact force with velocity, the vertical contact force declined significantly after the tire made contact with water flows. When the tire speed was 70 km/h, the rolling tire had not contacted water flows yet at 0.005 s, and the vertical contact force was equivalent to the applied loads. When the velocity reached 130 km/h, the aquaplane phenomenon occurred, and the vertical contact force was about 3000 N, which was about 25% lower than that under drying conditions. Furthermore, it can be found that along with the increase of tire speed, all the pavements showed a decreasing trend. OGFC could provide the tire with greater force on the vertical contact direction, followed by SMA and AC. In other words, the water-drainage properties and pavement texture could account for the different anti-skid performance of pavement. Additionally, as the contact duration reached 0.02 s, it also can be found that the tire had basically reached a balanced state of force. When the velocity increased from 70 km/h to 130 km/h, the ultimate tire–pavement vertical contact force decreased by 40~60%.

**Figure 20.** Effects of tire rolling velocity on vertical contact force. (**a**) AC pavement; (**b**) SMA pavement; (**c**) OGFC pavement.

#### *3.3. Effects of the Thickness of Water Film on Partial Aquaplane Performance*

With the other constant simulation parameters, the tire speed was set at 70 km/h for the investigation of the effect of water film thickness, that is, 4 mm, 8 mm, and 12 mm, on the vertical contact force under partial aquaplane circumstances. Based on historical output, the variation trends of tire–pavement vertical contact force on AC, SMA, and OGFC in partial aquaplane process with thickness of water film were obtained, as shown in Figure 21.

**Figure 21.** Effects of the thickness of water film on vertical contact force. (**a**) AC pavement; (**b**) SMA pavement; (**c**) OGFC pavement.

It can be seen from Figure 21 that while the thickness of the water film was 4 mm, the decreasing trends of vertical contact forces on the pavements were relatively stable. Furthermore, the ultimate vertical contact force decreased to about 50% of the applied loads onto the tire. With the increase of thickness of the water film, the decreasing trends of vertical contact force became obvious. As the thickness of water film reached 12 mm, the ultimate vertical contact force declined significantly, accounting for only about 15% of the applied vertical loads. In addition, while in the given thickness of water film, the vertical

contact force of OGFC was still superior to SMA and AC, which is still due to the difference of pavement textures.

#### **4. Conclusions and Future Research**

In this paper, steady state rolling analyses of the tire under free rolling state and braking state with a slip rate of 20% were conducted based on the finite element analyzing method in ABAQUS. Influences of pavement texture, tire pressure, and loads on tire–pavement friction behaviors were briefly discussed under drying conditions. Next, a tire–texture pavement contact model was built. The calculation results of steady state analysis were input into the explicit analysis by compiling INP files, and then the Euler grid model was constructed. Material parameters of water were defined, and the aquaplane finite element models of a radial tire on three types of asphalt pavement were constructed using the CEL method. Finally, variation laws of the tire–pavement vertical contact force with pavement type, tire rolling velocity, and thickness of water film under the partial aquaplane condition were analyzed by the built aquaplane model. Some conclusions were drawn as follows.

Based on the conclusions drawn above, FEA method, accompanied with laser scanning for the acquirement of pavement profile, seems to be a feasible method for analyzing tire– pavement friction on partial tire aquaplane conditions. Finally, a few future research topics are recommended as follows.

**Author Contributions:** Conceptualization, Z.Y.; Data curation, M.Y., Y.K. and L.Y.; Funding acquisition, M.Y. and L.Y.; Methodology, M.Y., J.L. and L.K.; Project administration, M.Y.; Software, Y.K. and J.L.; Writing—original draft, Y.K.; Writing—review & editing, Z.Y. and L.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is financially supported by Open Fund of Key Laboratory of Road Structure and Material of Transport Ministry, Chang'an University, the Fundamental Research Funds for the Central Universities, CHD (No.<PHONE_NUMBER>06), China Postdoctoral Science Foundation (No. 2018M633444), Shaanxi Province Postdoctoral Science Foundation (No. 2018BSHEDZZ123), General Project of Chongqing Natural Science Foundation (No. cstc2021jcyj-msxmX0554), National Natural Science Foundation of China (No. 51608085), and Key R & D Program of Guangxi (No. GUI KE AB20159036).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Characterization of Rutting Damage by the Change of Air-Void Characteristics in the Asphalt Mixture Based on Two-Dimensional Image Analysis**

**Kang Zhao 1, Hailu Yang 1, Wentao Wang 1 and Linbing Wang 2,\***

**Abstract:** In the process of the rutting test, the air-void characteristics in asphalt mixture specimens are a dynamic change process. It is of great significance to systematically study the correlation between the change of air-void characteristics and the depth of the rutting slab and establish a relationship with damage. In this paper, the air-void information of rutting specimen sections with different loading cycles (500, 1000, 1500, 2000, 2500, and 3000 times) is obtained by two-dimensional image technology. The dynamic change process of the micro characteristics of internal air voids of two graded asphalt mixtures (AC-13 and AC-16) under cyclic wheel load is analyzed, and it is used as an index to characterize the microstructure damage of the asphalt mixture. The results show that the variation of air-void distribution, air-void shape characteristics, and air-void fractal dimension with the loading process can well characterize the permanent deformation law of the rutting slab. The fractal dimension of the air void increases with the increase in load. It is a dynamic process in which the air-void content changes with crack initiation and propagation. After rutting deformation, the total air-void area and average air-void size of the sample increase, and the total air-void number decreases. Because microcracks are formed in the specimen after rutting damage, the aspect ratio of the air void increases, and the roundness value decreases.

**Keywords:** air-void characteristics; fractal dimension; crack initiation and propagation; rutting damage

#### **1. Introduction**

An asphalt mixture is a composite material composed of aggregates, air voids, mastic, and other multiphase media, which is generally considered a continuous medium. However, the anisotropy of the microstructure leads to complex mechanical states and complex interactions between the components of the asphalt mixture. Among these, the air-void characteristics significantly impact the mechanical properties of the asphalt mixture. We found that two specimens with the same gradation and same air-void content may have different air-void numbers and average air-void size distributions. Specimens with the same air-void content may be composed of air voids of different sizes according to a certain combination. Therefore, two specimens with the same air-void content may also have different failures under the same load conditions. Specimens with larger air voids are more prone to relatively early failure compared with specimens with smaller air-void sizes. This is because once the specimen is subjected to load conditions, the possibility of inducing larger strains increases. In turn, it may lead to internal structural instability, and the existing air voids may continue to expand in the specimen and connect (or merge with other air voids), resulting in serious damage [1]. The increase in air void will reduce the resistance of the asphalt mixture to pavement damage. This is because the air void cannot transfer the load, and the material becomes weaker due to the decrease in the effective area (The

**Citation:** Zhao, K.; Yang, H.; Wang, W.; Wang, L. Characterization of Rutting Damage by the Change of Air-Void Characteristics in the Asphalt Mixture Based on Two-Dimensional Image Analysis. *Materials* **2022**, *15*, 7190. https:// doi.org/10.3390/ma15207190

Academic Editor: Giovanni Polacco

Received: 24 August 2022 Accepted: 13 October 2022 Published: 15 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

concentrated force of each area is higher). The air-void content of the asphalt mixture shall be within a reasonable range to avoid adverse damage such as rutting [2]. Therefore, the proper compaction degree is significant for the performance of asphalt materials. The compaction index is widely used in engineering practice to ensure good service performance of asphalt mixture [3,4].

In recent years, with the progress of technology and the improvement of equipment available, the digital imaging technology of asphalt concrete is an effective tool for evaluating the internal structure. Researchers have been able to use digital image technology to monitor the internal failure process of engineering materials and connect this failure with the measured strain. Internal failure in materials can be represented in many forms, such as specific voids [5], crack surfaces [6,7], and spacing between cracks [8,9]. Air voids play an important role in characterizing the performance of asphalt mixture, and their distribution is very important for determining the overall mechanical response of asphalt mixture [10–12]. Under load, existing air voids may merge, resulting in microcracks at the interface between aggregate and mastic. Microcracks continue to expand and grow under deformation, forming macro cracks and increasing air voids [13]. Xu G et al. [14] extracted the internal void structure characteristics of asphalt mixture through an X-ray computed tomography (CT) test and three-dimensional (3D) image reconstruction technology. The change of air-void distribution before and after the freeze–thaw test is analyzed to evaluate the structural evolution of materials under freeze–thaw cycles. Yang B et al. [15] studied the correlation between the performance of porous asphalt mixture and the three pore characteristic parameters obtained by the CT. It was found that the micro pore characteristic parameters had little effect on high temperature, humidity sensitivity, and cooling performance but had a strong correlation with spalling resistance, permeability, connectivity, noise reduction, and other properties. Zhang Z et al. [16] used the connective void content of three types of asphalt mixtures was employed to characterize the damage to the corresponding asphalt mixture sample under freeze–thaw cyclic loading. The variation of connective void content revealed the nonlinear characteristics of asphalt mixture damage accumulation. Kassem Emad et al. [17] used several mechanical tests (the overlay tester, Hamburg Wheel-Tracking Test (HWTT), and a repeated tensile test) to characterize the influence of air-void distributions on mechanical properties and response of asphalt mixtures. The results show that air voids play an important role in influencing the performance of asphalt mixtures. Xu H et al. [18] used a set of image analysis programs to extract the internal structural characteristics of asphalt mixture during freeze–thaw cycles. The evolution of the internal void structure of the asphalt mixture during the freeze–thaw cycle was evaluated. Hassan et al. [19] analyzed the air void and crack characteristics caused by stress and strain by using a two-dimensional (2D) image analysis method and took them as the damage index to characterize the microstructure damage of the asphalt mixture. Therefore, it is very meaningful to analyze the damage behavior of asphalt mixture by describing the change characteristics of air voids in the deformation process of samples. However, previous researchers used image analysis to characterize the void of asphalt mixture, focusing on the change value before and after the test and less on the development of the whole process of void during the test. At the same time, the change of air-void content is more used as an evaluation indicator, and other intuitive and effective indicators are lacking. Few studies have used the rutting test to characterize the change in the internal voids of materials.

The purpose of this paper is to quantify the change of air-void characteristics as a damage characterization method by analyzing the section image information of rutting specimens during the test. This method can not only be used to quantify the changes of air-void characteristics in the process of damage development but also to determine the damage concentration area and use the determined air-void parameters to describe the severity of the damage. Although most of the current research involves 3D analysis, the internal structure of materials is obtained through various advanced scanning techniques. It is also important to understand the basic 2D measurement in damage analysis because it

can provide a reference for specific problems. Using the same concept, further analysis of damage can be extended to complex problems in 3D or four-dimensional (4D) analysis [20]. At present, most of the research focuses on the change of meso-void characteristics of the whole specimen before and after wheel loading. Few studies have explained the airvoid distribution of the whole wheel load area section in the whole process of loading and linked the air-void distribution with the change of air-void fractal dimension and permanent deformation of the asphalt mixture. This study is a comprehensive study on the variation law of internal air-void characteristics of two kinds of asphalt mixtures (AC-13 and AC-16) with the loading cycles of the rutting test. It is the basic work to explain the evolution of internal defects of mixtures under wheel load and evaluate the high-temperature performance of asphalt mixtures from a meso perspective.

#### **2. Materials and Methods**

#### *2.1. Materials*

In this study, the dense gradation, namely AC-13 and AC-16, with the nominal maximum aggregate size of 13 and 16 mm, was selected. Qinhuangdao 70# asphalt was selected as an asphalt binder, and the main physical properties are given in Table 1. The coarse aggregate is limestone, the fine aggregate is machine-made sand (0–5 mm), and the filler is limestone powder. Asphalt mixture samples are prepared under the optimum asphalt binder content obtained by the Marshall method [21]. The aggregate gradation, optimum asphalt binder content, and air-void content of the two asphalt mixtures are shown in Table 2. The air-void content of all types of samples is controlled within ±0.5% of the range shown in Table 2.

**Table 1.** The properties of asphalt binders.

**Table 2.** Aggregate gradations and mix design results of AC mixtures.

#### *2.2. Specimen Fabrication and Rutting Test*

According to the Chinese code [22], the mixing temperature of asphalt mixture containing petroleum asphalt is 163 ◦C. The compaction temperature of the rutting specimen is 100 ◦C during the forming process. Therefore, the rolling equipment and test mold shall be preheated to 100 ◦C in advance. The roller load is set as 9 kN (line load 300 N/m). The manufacture of a rutting specimen consists of two steps: (1) rolling the rutting specimen for two cycles with a roller (moving back and forth on the sample is defined as one cycle) and (2) rolling the sample for 12 cycles after rotating it 180◦. The size of the final formed rutting specimen is 300 (length) × 300 (width) × 50 mm (height). In order to accurately obtain the change of specimen section after different loading cycles, the plate is cut into three parts by an asphalt mixture cutter along the wheel loading direction before the test, 60 mm in the middle and 120 mm on both sides. Therefore, 3 specimens named 1, 2, and 3 are obtained, as shown in Figure 1.

**Figure 1.** Specimen preparation and cutting.

The rutting test is carried out by using the automatic asphalt mixture rutting instrument to study the rutting development of the asphalt mixture. The device is mainly composed of rubber wheels and a fixed platform. The diameter of the rubber wheel is designed to be 200 mm, and the contact width of the wheel specimen is 50 mm. The stress applied by the rubber wheel is 0.70 ± 0.05 MPa, and it is loaded back and forth on the rutting specimen at the rolling speed of 42 r/min. The experimental temperature is set at 60 ◦C [22]. Due to the cutting effect, the overall size of the rutting specimen is reduced by 4 mm. In order to keep the boundary conditions unchanged, 4 mm steel sheets are inserted on one side to provide a constant lateral force. The details of the rut test are shown in Figure 2. The cumulative rutting depth changes of the specimens under different loading times (500, 1000, 1500, 2000, 2500, and 3000) are recorded by the linear variable differential transformer (LVDT) in the device. In this study, 3 rut specimens of AC-13 and AC-16 are used, respectively.

**Figure 2.** Rutting test.

#### *2.3. Image Acquisition and Processing*

The industrial camera is used to collect the images of four sections of each specimen after the corresponding loading times (500, 100, 1500, 2000, 2500, and 3000) to study the rutting deformation law. Because of the large difference in each specimen section, we only collected the image information of one of the three samples of each asphalt mixture. When taking images, the camera is placed on the test bench to ensure the levelness of the camera. The camera lens is located at the center of the specimen, 50 cm away from the specimen, as shown in Figure 3.

**Figure 3.** Specimen section image acquisition process.

Digital image processing uses a series of algorithms to process and analyze digital images with computers so that the images can meet the needs of human vision, other equipment, data extraction, etc. [23]. In this section, MATLAB software is used to process images, including image enhancement, image filtering, image segmentation, and feature information acquisition. (1) Image enhancement: image enhancement is mainly processed by histogram equalization, also known as gray-level equalization. The purpose is to convert the input image into an output image with the same number of pixels at each gray level through the operation of each pixel in the image. The cross-section image of the asphalt mixture specimen after gray histogram equalization is shown in Figure 4b. Histogram equalization function in MATLAB is histeq (f, 256). (2) Image filtering: In this study, median filtering method is used to filter the image, which is very effective for filtering pulse interference and image scanning noise. The cross-section of the specimens after median filtering is shown in Figure 4c. The function of the median filtering operation in MATLAB is medfilt2 (f). (3) Image segmentation: This research obtains the best threshold value through the function gray threshold value in MATLAB [24]. According to the function result, when the threshold value is set to 35, a good segmentation effect can be obtained, as shown in Figure 4d. Figure 4e shows the acquisition of air-void feature information on the segmented image using ImageJ software. The whole process is shown in Figure 4. By comparing with the air-void content of the actual mixture, we selected the air void when the intensity of the corresponding pixel is in the range of 0 to 20. The air voids are divided into small air voids (<1 mm2), medium air voids (1–5 mm2), and large air voids (>5 mm2) according to the size.

#### *2.4. Air-Void Characteristics*

Some parameters, such as air-void content, air-void number, air-void shape, and average air-void size, have been successfully used to characterize the characteristics of air voids. Among them, the number and average size of air voids are very important in characterizing the characteristics of air voids. In the study of microstructure damage, the changes of these parameters (by comparing the parameters during and before deformation) provide information on the damage degree of rutting specimens with different gradations, different depths, and different loading accumulation times.

**Figure 4.** Image processing process: (**a**) original image; (**b**) section image after gray histogram equalization; (**c**) section image after median filtering; (**d**) sectional image after image segmentation; (**e**) air-void feature information.

In 2D images, the shape is usually regarded as the area surrounded by a closed contour curve. The most commonly used shape factors in image analysis are roundness, circularity, and aspect ratio. Figure 5 shows air-void shape features defined according to roundness, circularity, and aspect ratio. The shape factor is usually normalized, ranging from 0 to 1. A shape factor equal to 1 usually represents the ideal case or maximum symmetry, such as a circle and sphere [25]. Generally speaking, the crack has a high aspect ratio, and the circularity and roundness values are low (close to 0.0), showing a flat shape. In this study, the section information of the specimen before and after the rutting test is obtained through image analysis. The average air-void size, aspect ratio, circularity, and roundness of the air void in the whole process of the test are calculated by using Equations (1)–(4) to study the change of the air void.

$$\text{Average air voids size} = \frac{\text{Total air void area}}{\text{Number of air void}'} \tag{1}$$

$$\text{Asspect ratio} = \frac{\text{Length of major axis}}{\text{Length of minor axis}'} \tag{2}$$

$$\text{Circularity} = \frac{4\pi \times \text{Area}}{\text{Perimeter}^2} \text{\textdegree} \tag{3}$$

**Figure 5.** Shape features.

#### *2.5. Fractal Dimension Theory*

The fractal dimension reflects the validity of space occupied by complex objects. It is a measure of the complex shape and irregularity of objects, including Hausdorff dimension, box-counting method, etc. [26]. Box counting uses a set of square boxes or grids to measure the length or distance between points on the shape boundary. It can be calculated with MATLAB software [27]. Figure 6 shows an image of a box array of different sizes (r) covering the air-void image. In the curves of Log (box number, N) and Log (box size, r), the boxes containing the pixels of the air-void image are counted (N), which is used to obtain the fractal dimension, i.e., the slope of the logarithmic regression line. As shown in Figure 7, it can be seen from the fitting curve of Log N (s)–Log N that the image shows the self-similar characteristics of the fractal body at different scales. This shows that the fractal dimension based on 2D images can well reflect the changes in micro-structure. Therefore, it is feasible to quantitatively characterize the rut section by fractal dimension. We can use Equation (5) to calculate the fractal dimension. Lower fractal dimensions represent more stable air-void structures. If the measured air-void areas are the same, different air-void distributions will produce different fractal dimensions.

$$FD = -\frac{\log(N(r))}{\log(r)},\tag{5}$$

**Figure 6.** Images of box arrays of different sizes (r).

**Figure 7.** Logarithmic relationship between box count and each box size.

#### **3. Results**

#### *3.1. Correlation between Rutting Performance and Air-Void Characteristics*

#### 3.1.1. Rutting Test Results

Rutting depth is an important index to evaluate the rutting resistance of specimens. The average value of rutting depth measured on different specimens (three samples for each grading) of the same grading is taken as the rutting test results of this study. The rutting test results are shown in Figure 8. As shown in Figure 8, the deformation of AC-16 is less than that of AC-13 during the whole flow deformation. Different asphalt mixture gradation leads to different deformation development of rutting specimens. Later, we try to explain this phenomenon from the damage analysis of the screenshot image of the rutting specimen.

**Figure 8.** Rutting test results.

#### 3.1.2. Correlation between Rutting Cumulative Depth and Air-Void Characteristics

Previous studies have found that the damage of the asphalt mixture under load may first occur in its internal air voids [28]. Therefore, it is of great significance to deeply understand the permanent deformation of the specimen by studying the variation characteristics of the air voids in the specimen during the rutting test. The correlation between the accumulated rut depth and the air-void characteristics (average air-void size, air-void content, circularity, and aspect ratio) is shown in Figure 9. In order not to increase the influence of

other factors, only the air-void characteristics of section 2-A after different loading times are studied here.

**Figure 9.** Correlation between rutting depth and air-void distribution: (**a**) average air-void size; (**b**) air-void content; (**c**) circularity; (**d**) aspect ratio.

The cumulative rut depth is linearly related to the average air-void size and aspect ratio, as shown in Figure 9a,d. With the increase in the average air-void size and aspect ratio, the cumulative rut depth also increases. At the same time, according to Figure 9c, with continuous loading, the air void inside the specimen becomes flatter, and the circularity decreases. At this time, the internal damage of the specimen increases, and it is more prone to damage. It can be seen from Figure 9b that with the increase in cumulative rut depth, the air-void content first increases, then decreases, and finally tends to a stable value. This is mainly because the number of air voids will first increase rapidly under the load. With the increase in load, the air voids in the specimen are continuously compacted and merged into larger air voids, finally forming cracks. Therefore, the air-void content decreases to a certain extent. This is consistent with Ma X et al.'s research on two possible air-void evolution laws in the mixture [29]. This indicates that under load, the change of internal air voids, along with the initiation and propagation of cracks, is a dynamic process, which is of great significance to the study of internal damage of rut specimens.

#### 3.1.3. Correlation between Rutting Cumulative Depth and Air-Void Fractal Dimension

Based on the fractal theory, the relationship between the fractal dimension of air voids and the rutting depth of rutting specimens is studied. It solves the difficulty of measuring the deformation of each layer in the sample by a conventional test method. The measured fractal dimensions of the AC-13 and AC-16 sections vary with the loading cycles, as shown in Figure 10. For the relationship between fractal dimension and rutting depth, we only analyze the fractal results of section 2-A, and the results are shown in Figure 11. The fractal dimension is related to the complexity and disorder of surface morphology, and more complex image surfaces have higher fractal dimensions [30]. In general, the fractal dimension of air voids increases with the increase in loading times. In other words, with the increase in loading times, the air-void complexity increases. This is consistent with the formation of many irregularly shaped air voids and the increase in roughness in the specimen section during the test.

**Figure 10.** Relationship between fractal dimension and loading cycles: (**a**) AC-13; (**b**) AC-16.

**Figure 11.** Relationship between fractal dimension and rutting depth.

By studying the correlation between rutting depth and air-void fractal dimension, it is found that the value of the fractal dimension is linearly correlated with rutting depth. Therefore, fractal theory can reliably study the deformation law of rutting samples. With the development of deformation, the fractal dimension of cross-section air voids increases gradually. During the test, air voids of all sizes are compacted under the wheel load and create new air voids, especially in the area below the load. With the development of rutting depth, new air voids appear, and the old air voids fuse and merge into micro cracks, which increases the complexity and dispersion in the section. Obviously, the fractal results of air voids accord with the physical significance of fractal theory.

#### *3.2. Air-Void Structure and Rutting Specimen Damage*

#### 3.2.1. Change of Air-Void Content

In order to characterize the internal damage of the material by the change of the air-void content of the specimen. The changes in the air-void content of the 2-A section of the two mixtures with the loading cycles and the height along the specimen before and after loading are studied, respectively, as shown in Figure 12.

**Figure 12.** Change of air-void content: (**a**) change with loading cycles; (**b**) along the depth direction.

The air-void content of the two rutting specimens first increases, then decreases, and finally stabilizes with the increase in loading cycles, as shown in Figure 12a. This is mainly due to the formation of a large number of air voids in the loading area under compression and shear at the early stage of loading. The newly formed air voids are more than the air voids used for fusion and forming microcracks, so the air-void content will have a rising process. With the increase in loading cycles, a part of the internal air void of the rutting specimen is fused and merged to form micro cracks under load, and a part of the air void is compacted, so the internal air-void content will have a decline process. In the last stage, damage occurred inside the specimen. The change process of air voids is similar to the previous two stages, but this process is not as violent as before. Finally, the air voids will gradually tend to a stable value. Due to the lack of laboratory technology to accurately capture the air-void distribution of asphalt mixture, previous research results are based on the air-void content before and after the test [31]. It is of little significance to study the damage development process.

It can be seen from Figure 12b that before and after the test, the area with high air-void content change mainly occurs in the middle and upper part, which indicates that under the action of wheel load, this area is subject to more shear deformation, resulting in more new air voids. For AC-13, the air-void content in the middle area of the section decreases after the test. Because this area is mainly subjected to compressive stress and the air voids are continuously compacted, the air voids are reduced. The increase in bottom air-void content indicates that it is subjected to shear stress. For AC-16, the change of air-void fraction in the middle area is not obvious, and there is a certain decrease in the bottom. This is mainly due to the large particle size of aggregate in AC-16 and the good performance of the skeleton in the middle area, which can withstand large shear stress failure.

#### 3.2.2. Damage Analysis of Rutting Specimen

Previous studies have shown that the load transfer capacity of asphalt mixture mainly depends on the interaction between aggregate and mortar, and macro cracks appear at the interface between aggregate and mortar [32]. Figure 13 shows the change process of airvoid shape of the two mixtures under wheel load and also vividly represents the formation

process of macro cracks in the mixture. As shown in Figure 13b, under the action of wheel load, the air void develops along the interface between aggregate and asphalt mortar. This process leads to the deformation, expansion, and fusion of the air void and promotes the generation of macro cracks in the mixture. Micro damage also promotes the bond failure or damage of asphalt mortar, further has a more adverse impact on the bonding performance between aggregate and asphalt mortar, and aggravates the degradation of internal structure under a wheel load. In addition, it can be seen from Figure 13a that the load also makes the aggregate with smaller particle size in the mixture crushed, which destroys the internal skeleton structure of the mixture and further accelerates the damage of the asphalt mixture. Comparing the two mixtures, it can be seen that this phenomenon is more obvious for the mixture with a smaller maximum particle size of aggregate. The damage is not only the increase in air void but also accompanied by the crushing of many aggregates and the overflow of asphalt mortar. This is different from Li P et al.'s result that only aggregate movement exists in asphalt mixture under pressure [33], which provides a new idea for the study of mixture damage.

**Figure 13.** Variation of sectional air void with loading cycles: (**a**) AC-13; (**b**) AC-16.

3.2.3. Change of Air-Void Characteristics

Figure 14 shows the changes in total air-void number, total air-void area, and average air-void size of the two graded mixtures before and after the rutting test. These parameters have been successfully applied to characterize the characteristic changes of voids and have been well-established in previous research work [34,35]. In the study of microstructure damage, the changes of these parameters (by comparing the parameters before and after deformation) provide valuable information on the severity of damage at different depths in compacted asphalt mixture specimens. The bar graph compares changes in air-void characteristics by averaging at 5 mm intervals along the height of the specimen. Through observation, several noteworthy results were obtained along the characteristic distribution of the sample height in three main areas, namely, the bottom (0–1.5 cm), the middle (1.5–3.5 cm), and the top (3.5–5 cm).

**Figure 14.** Variation of air-void characteristics in rutting test: (**a**) AC-13; (**b**) AC-16.

As shown in Figure 14, the rate of change of total air-void area and average air-void size is higher in the middle and upper region, which is consistent with the slight expansion observed in the middle and upper region of the specimen after the rutting test. Under the action of wheel load, the middle and upper mixture is seriously damaged. The total air voids of the two mixture specimens along the height direction are reduced, indicating the generation of micro cracks and macro cracks in the rutting specimen, which also promotes the expansion in the middle area of the specimen. This is mainly due to two reasons: one is that the internal air void is reduced due to the continuous compaction of the specimen under load; Second, microcracks are formed due to the fusion of small air voids.

In this study, when describing the different damage phenomena in the specimen due to deformation, the changes in total air-void number and total air-void area are determined as two key parameters. According to the previous study [36], there are two main mechanisms under permanent deformation, namely densification, and shear deformation, as shown in Figure 15. By analyzing the changes in air-void number and air-void area before and after the test, the results of this experiment are mainly related to the damage mechanism related to compression deformation. The relationship between damage mechanisms and air-void properties will be further clarified in subsequent sections, involving changes in air-void geometry.

**Figure 15.** Correlation between damage mechanisms and changes in air-void content properties.

The change distribution of air-void characteristics along the height of the rutting specimen before and after the test is systematically studied in the previous part. This part will study the change of air void under load from the whole section. Firstly, the variation characteristics of air-void perimeter and area in the mixture under different loading cycles are studied. Secondly, the changes in air voids with different area sizes and their corresponding quantities before and after loading are studied. Finally, the changes in the total area and number of air voids of different types (large, medium, and small) before and after the test are studied.

Figure 16 shows the area and perimeter distribution of the air void in the rutting specimen under different loading cycles. When the polygon area is certain, the larger the ratio of major and minor axes, the longer the circumference [37]. It can be seen from Figure 16 that with the increase in loading times, the area and perimeter distribution of the internal air voids of the rutting specimen become wider, and the discreteness also increases. This shows that under the action of load, the air voids develop towards a flatter shape, preparing for further integration and macro cracks. This also confirms the rationality of the previous research on the development law of air voids before and after loading from another perspective.

**Figure 16.** Relationship between air-void area and perimeter: (**a**) AC-13; (**b**) AC-16.

Figure 17 shows the variation law of the average area and perimeter of air voids along the height of the specimen. It can be seen from Figure 17 that before the load is applied, the average perimeter and area of the two graded mixtures change a little along the height direction, which shows that the specimen is fully compacted during the preparation process. After 3000 cycles of cyclic loading, the average area and perimeter of the specimen increase to a certain extent along the height direction, indicating that the internal air-void size of the specimen increases as a whole. For the two mixtures, the change rate of the middle and upper part of the specimen is large, which is also consistent with the obvious expansion of the middle and upper part observed after the test.

In order to further explain the change of air voids before and after the test, we calculated the area and the corresponding number of air voids in the 2-A section of the specimen before and after loading, as shown in Figure 18. Previous field tests have concluded that air voids in the pavement are compressed after permanent deformation [38]. However, it can be seen from Figure 18 that a larger area of air voids appeared in the specimen after loading. The number of air voids with smaller areas (≤1 mm2) is reduced greatly, and most of them are compressed. At the same time, we also compared the change rate of the total area and quantity of different types of air voids before and after the test, as shown in Figure 19.

**Figure 17.** Variation of the average area and perimeter of air voids along specimen height (**a**) Average area; (**b**) Average perimeter.

**Figure 18.** Relationship between air voids and corresponding quantities: (**a**) AC-13; (**b**) AC-16.

**Figure 19.** Changes in area and quantity of air voids with different sizes.

It can be seen from Figure 19 that after 3000 cycles of wheel loading, the total number and area of small air voids in the specimen decrease, while the middle air voids and large air voids increase, and the change rate of the total number and area of large air voids is the largest. This is mainly due to the small number of large air voids in the original specimen. After loading, the number of large air voids increases, making the change rate larger. The above results show that the air voids inside the specimen mainly occur through compaction and fusion of small air voids under load. The fusion between small air voids expands into micro cracks, which in turn induces damage to the asphalt mixture. At the same time, the change rate of the total number and area of the large and medium air voids in AC-13 specimens is higher than that of AC-16, indicating that the damage degree of AC-13 specimens is higher than that of AC-16. This result can be confirmed by rutting results.

#### 3.2.4. Change of Air-Void Shape Characteristics

After the rutting test, not only the total number of air voids, total air-void area, and average air-void size but also the shape characteristics of air voids will change greatly. Figure 20 shows the change rate of void shape characteristics of the specimen before and after the test, i.e., aspect ratio, circularity, and roundness. The bar graph compares the changes in air-void shape characteristics by averaging at 5 mm intervals along with the height of the specimen.

**Figure 20.** Variation of air-void shape characteristics in rutting test: (**a**) AC-13; (**b**) AC-16.

It can be seen from Figure 20 that after the test, the overall roundness of the air voids in the specimen decreased, and the shape of the air voids developed into a flat shape. At the same time, with the increase in air-void size, two adjacent air voids merge to form an air void, and its aspect ratio is higher than that of its single air void. This confirms the previous findings about the change in air-void properties. In general, the method of measuring the shape characteristics of air voids provides a different method to describe the characteristics of air-void formation, growth, connectivity, and expansion. The damage is due to the coalescence of small air voids to form microcracks, which increases the aspect ratio of air voids, thereby reducing the number of air voids and roundness [39]. Based on the observation of a large number of samples, it is determined that the damage controlled by the change of air-void structure occurs in many ways. First of all, these changes may be due to the increase in the existing single air-void size and the propagation of cracks. Secondly, the two separated air voids may merge and expand, forming cracks near the damage. Third, due to deformation, new air voids with crack initiation potential are also formed. This finding is consistent with the two failure types determined by Kim et al. [40]. The damage in the asphalt mixture is the result of microcracks generated in the interfacial transition (ITZ) between asphalt and aggregate and within the binder (cohesive failure).

#### **4. Conclusions**

This paper aims to quantify the change of air-void characteristics as a method to characterize damage by analyzing the cross-sectional image information of rutting specimens during the test. Firstly, the indoor rutting test was carried out, and the air-void information of the sample interface under different loading times was obtained by digital image technology. Then, the variation law of the internal air-void characteristics of the asphalt mixture with the loading times of the rutting test is systematically studied, and the relationship with the development of rutting is established. Finally, the air-void characteristic parameters are used to quantify the damage in the asphalt mixture and describe the generation and development of damage. According to the research results, the following conclusions can be drawn:

Therefore, we can use the change of void ratio, void characteristics, and void shape characteristics to quantify the damage type and its evolution process in the rutting formation of the asphalt mixture in future research. Due to the limitations of the technology used in this study, the voids discussed are obtained based on two-dimensional images. It is difficult to distinguish the connection and closure of voids in the three-dimensional case of a mixture, so further research is needed. In addition, in order to determine the main cause of damage to the inner cup of the asphalt mixture, it is necessary to analyze the damage at the contact between aggregate and asphalt and the impact of the contact point between aggregates on the performance.

**Author Contributions:** Conceptualization, K.Z. and L.W.; methodology, L.W. and W.W.; software, K.Z.; resources, H.Y.; data curation, K.Z.; writing—original draft preparation, K.Z.; writing—review and editing, L.W.; supervision, L.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Beijing Major Science and Technology Projects, grant number (No. Z191100008019002).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

#### **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Study on the Extraction of CT Images with Non-Uniform Illumination for the Microstructure of Asphalt Mixture**

**Lei Zhang 1, Guiping Zheng 1, Kai Zhang 2, Yongfeng Wang 3, Changming Chen 3, Liting Zhao 3, Jiquan Xu 4, Xinqing Liu 4, Liqing Wang 4, Yiqiu Tan 1 and Chao Xing 1,\***

**Abstract:** An adaptive image-processing method for CT images of asphalt mixture is proposed in this paper. Different methods are compared according to the error analysis calculated between the real gradation and 3D reconstruction gradation. As revealed by the test results, the adaptive imageprocessing method was effective in carrying out different brightness homogenization processes for each image. The Wiener filter with 7 × 7 size filter was able to produce a better noise reduction effect without compromising image sharpness. Among the three methods, the adaptive image-processing method performed best in the accuracy of coarse aggregate recognition, followed by the ring division method and the global threshold segmentation method. The error of the gradation extracted by the adaptive image-processing method was found to be lowest compared with the real gradation. For a variety of engineering applications, the developed method helps to improve the analysis of CT images of asphalt mixtures.

**Keywords:** asphalt mixture; Computed Tomography; image process; microstructure; gradation analysis

### **1. Introduction**

As variable multi-phase composite material, asphalt mixture is characterized by significant differences, randomness, and variability among the constituent materials, which contributes to the non-uniqueness of its microstructure. However, the spatial distribution characteristics of each component in the internal structure of asphalt mixture cannot be reflected by the traditional macro-empirical evaluation method, which focuses primarily on the overall features and road performance of asphalt mixture at the macro scale. However, for multi-phase asphalt mixture, the macro mechanical properties depend on its microstructures. In this circumstance, it is necessary to calculate the performance mechanism of asphalt mixture from a meso-level perspective. At present, computed tomography (CT) and digital image-processing technology have been widely applied in micro-scale research on asphalt mixture. It is difficult to properly identify the differential materials in CT images of asphalt mixture, due to factors such as CT imaging mechanism and variations of asphalt-mixture gradation composition. Meanwhile, grayscale values of substances with the same density should theoretically be consistent at different locations, but the scanning energy limitation of industrial CT machines causes them to present varying degrees of lightness and darkness at different distances from the specimen's central axis. Therefore, it remains challenging to accurately identify aggregates, voids, and asphalt mortar. To obtain an accurate extraction of the microstructure of asphalt mixture, it is necessary to further study the image-processing method used for CT images of asphalt mixture.

It was in the early 1990s that research began on the CT image-processing method used for asphalt mixture. At that time, Masad et al. [1] applied CT technology and digital camera technology for the first time to conduct a quantitative analysis of the compaction effect of

**Citation:** Zhang, L.; Zheng, G.; Zhang, K.; Wang, Y.; Chen, C.; Zhao, L.; Xu, J.; Liu, X.; Wang, L.; Tan, Y.; et al. Study on the Extraction of CT Images with Non-Uniform Illumination for the Microstructure of Asphalt Mixture. *Materials* **2022**, *15*, 7364. https://doi.org/10.3390/ ma15207364

Academic Editor: Simon Hesp

Received: 5 September 2022 Accepted: 10 October 2022 Published: 20 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the mixture under different compaction modes, which led to success in introducing CT into the research on the microstructure of asphalt mixture, thus opening up a new field of research. In the meantime, research was advanced on the image-processing method applied for composite CT.

CT scanning technology is used to detect the attenuation of X-rays after passing through an object. Then, the attenuation images are outputted to the reconstruction algorithm for completing the reconstruction [2]. Due to the poly-chromatic nature of X-ray energy, however, the attenuation coefficient can cause errors, thus resulting in artifacts within the reconstructed images [3,4]. The methods used to correct artifacts can be categorized into physical filtering methods and software correction methods. Commonly used during scanning, the physical filtering method [5–8] can reduce the artifacts. Meanwhile, however, it also reduces ray intensity, which contributes to SNR (signal noise ratio) reduction. Due to fewer restrictions during application, software correction methods [3,9–11] are also frequently used. For asphalt mixture, however, which is a multiphase material, these methods are more complicated for correcting artifacts due to the lack of prior knowledge about X-ray spectra and material properties.

In recent years, CT scanning technology has been increasingly used for the microstructural characterization of asphalt [12], bituminous mortar [13], and asphalt mixture. Within research on asphalt-mixture CT image processing, the focus has mainly been placed on the methods emloyed to separate aggregate, mortar, and void within the image [14,15]. Currently, commonly used methods include the artificial threshold method [16–18], fuzzy C-means algorithm [19,20], Gaussian mixture model [21], and Otsu algorithm [22–25]. The most widely used is the Otsu algorithm. When Otsu is applied to images with artifacts, however, the segmentation effect is less than satisfactory. To solve this problem, Liu put forward the ring division method [26], which achieved higher segmentation accuracy. Nevertheless, it failed to determine accurately the size of the test piece and its position in the image. Additionally, the method was restricted to reducing the impact of artifacts on image segmentation when faced with the artifacts distributed along the radial non-linear section of the test piece, meaning it was unable completely to address the problem of artifacts.

In order to address the current shortcomings of image-processing methods, this paper proposes a brightening program for asphalt-mixture CT images that can help effectively eliminate the adverse effect of artifacts on image threshold segmentation. The processing effect based on this method was compared with that produced by the ring division method and the global threshold segmentation method, respectively. The research findings of this paper can contribute to more precisely identifying the composition and different materials of asphalt mixtures, accurately ascertaining the microstructure of asphalt mixture, and accomplishing a quantitative analysis of void distribution characteristics, aggregate homogeneity, and microscopic damage characteristics of asphalt mixture. This can provide a thorough understanding of the internal microstructure and properties of asphalt mixture, and improve the relevant test methods and protocols to further improve the performance of pavements. In addition, the method can also be applied to cement concrete, composite materials, and metal materials, so it has promising application potential.

#### **2. Materials and Methods**

#### *2.1. Materials*

In this study, andesite produced in Heilongjiang Province and limestone powder produced in Jilin Province were used. The relevant tests were conducted on stone materials, in line with the standards JTGE42-2005 [27] and JTGF40-2004 [28]. According to the test results, all indicators met the requirements. SBS polymer modified asphalt was used in this study, and the properties are shown in Table 1. The lignin fiber content of the SMA asphalt mixture was 0.3%.

**Table 1.** Properties of SBS modified asphalt.

#### *2.2. Gradation Design*

AC-16, OGFC-16, and SMA-16 test pieces were designed and fabricated according to the standard JTGF40-2004, with three replicates prepared for each sample. The gradation is shown below in Figure 1.

**Figure 1.** Gradation curves of three kinds of asphalt mixtures.

#### *2.3. Optimum Asphalt Content and Test Pieces Preparation*

The optimum oil content of the asphalt mixture was determined using the Marshall asphalt mixture design method. The optimum asphalt content of AC-16, SMA-16, and OGFC-16 was 4.5%, 6%, and 3.1%, respectively. The mixtures were compacted 120 times using the rotary compaction method. To obtain more pictures of the test pieces, the height of the test pieces after molding was set to about 150 mm. The volume parameters of the asphalt mixtures are shown in Table 2.

**Table 2.** Volume parameters of asphalt mixtures.

#### *2.4. Scan by Industrial CT*

The CT scanning equipment used for Marshall test pieces was a Phoenix micro focus industrial CT purchased by the School of Transportation Science and Engineering affiliated to Harbin Institute of Technology. The cone-beam scanning industrial CT system and the cone-beam filter back-projection reconstruction algorithm were both applied to reconstruct the faulty images. The scanning voltage/current was 195 kv/95 μA. The critical performance indicators are shown in Table 3. In this study, 150 mm corresponded to 1500 pixels in the image, so the resolution of the asphalt-mixture image recognition was 0.1 mm.

For the original model subjected to CT scanning and 3D reconstruction, the crosssectional images were exported at equal spacing, with the spacing value set to 0.15 mm, and each test piece was capable of exporting a total of 1000 images.

#### **3. CT Image Feature Analysis**

Taking SMA-16 as an example, there were roughly five peaks shown in the gray histogram of the asphalt-mixture CT image. Under normal conditions, each peak is represented in ascending order of gray scale: background, voids and artifacts, asphalt mortar, and aggregates. In practice, however, the artifacts appeared in the center of the test piece within the image, i.e., the brightness near the center was lower than near the edge, and the brightness along the radius of the test piece was nonlinear, which is largely attributed to ray hardening. Consequently, the gray value of some aggregate shifted to the wave peak of mortar. In case of artifacts being left unprocessed, the aggregate was extracted directly using the global threshold segmentation, causing the aggregate within the artifacts to be recognized as background removal, thus resulting in the loss of aggregate information as shown in Figure 2. This is where the difficulty lies when processing asphalt-mixture CT images.

**Figure 2.** Characteristics of CT images. (**a**) Original image; (**b**) gray histogram; (**c**) global segmentation.

#### **4. Image-Processing Method for Asphalt Mixture**

*4.1. Adaptive CT Image-Processing Method*

4.1.1. Obtain the Center and Radius of the Test Piece

It is important to determine the center and radius of the test piece for subsequent image-processing operations. During the actual CT scanning process, however, there is no way to guarantee that the test piece is located at the center of the scanning platform. Thus, the test piece can show different center positions in the images. Therefore, it is necessary to identify the center and radius of the test piece in the image.

Firstly, the matrix of the entire image was extracted into a rectangular coordinate system, with the row and column numbers of the pixel points representing their Y and X coordinates respectively.

Then, the image was segmented using a single threshold. A range of values was set for the radius of the circle surface is, starting from the minimum value, and a circular surface with that radius value was drawn and allowed to move freely over the binarized image. As shown in Figure 3, the maximum number of white pixels it can cover was recorded, and then the radius value was increased by one pixel length. The above steps were repeated until the radius of the circular surface reached the maximum value, or the maximum number of white pixels covered after increasing the radius value by one pixel length remained unchanged. At that time, the iteration process was discontinued.

**Figure 3.** Diagram of center and radius's determination. (**a**) Circular surface covering white pixels; (**b**) determined radius of the test piece; (**c**) determined center of the test piece.

Initially, following the method described above, the upper limit on the radius of the circular surface was excessively small. The upper limit was increased and iteration continued. In the second case, the maximum number of white pixels that can be covered by the circular surface with the value of adjacent radius was identical, suggesting that the smaller value was the radius of the test piece. Then, the circular surface was drawn with the radius of the test piece, so that it could be moved within the image. When it covered the maximum number of white pixels, the center of the circular surface overlapped with the center of the test piece, meaning that the center of the test piece could be determined.

Finally, the radius and center of test piece were used for segmenting the image to remove the redundant information while improving the processing efficiency.

4.1.2. Image Brightness Homogenization

Since the distribution of the artifacts was a circular range concentric with the test piece, the brightness adjustment operation was conducted by calculating the distribution of the grayscale of the test piece along the radius direction, with the procedures detailed as follows.

(1) Zero processing of the grayscale of voids. Since the extracted object used in this study was aggregated, it was necessary for mortar and voids to be removed as background. However, mortar cannot be removed due to the presence of artifacts, resulting in the grayscale of mortar and making the aggregates overlap each other. In order to prevent this from impacting the calculation of the statistical distribution of grayscale along the radius, multi-threshold segmentation and image subtraction [22] were first performed to remove the voids, as shown in Figure 4.

(**a**) (**b**)

(**c**)

**Figure 4.** Zero processing of grayscale voids. (**a**) Cropped image; (**b**) extraction of voids (white part); (**c**) voids' gray value return to zero.

(2) Calculating the distribution of grayscale along the radial direction. Based on the previously determined center of the test piece, the average gray values of the pixels on the approximate ring concentric with the test piece were calculated one by one. The expression of each approximate ring is shown in (1), while the approximate ring is shown in Figure 5a.

$$i - 1 < \sqrt{(x - x\_0) + (y - y\_0)} \le i \tag{1}$$

where *x*0 and *y*0 represent the horizontal and vertical coordinates of the center of the test piece, respectively; *x* and *y* indicate the horizontal and vertical coordinates of the pixel points on the *i*-th approximate ring, respectively; *i* denotes the outer diameter of the approximate ring, and the value ranges from one pixel length to radius.

**Figure 5.** Obtained gray distribution along the radial direction. (**a**) Approximate ring concentric with the test piece; (**b**) radial distribution of grayscale.

After determining the approximate ring that different pixels belong to, the average gray value of the pixels on each approximate ring was calculated, as shown in Formula (2).

$$Average\_i = Sum\_i / Num\_i \tag{2}$$

where *Sumi* represents the sum of the pixels' gray value with the gray value of [0, 180] on the *i*-th approximation ring; *Numi* indicates the number of pixels with the gray value of [0, 180] on the *i*-th approximation ring; and *Averagei* denotes the average gray value of pixels with the gray value of [0, 180] on the *i*-th approximation ring.

If the radius of the test piece is 503 pixels in length, 503 average values can be calculated. To avoid affecting the brightening operation, inclusion was limited to pixels with gray values greater than 0 but less than 180 before the average value was calculated. This is because the aggregates with large gray value and voids did not need to be increased in brightness, so they were excluded from the calculation. Figure 5b shows the radial distribution of gray value.

(3) Based on the average gray value of the pixels in each approximate ring, 180 was taken as the target gray value. When the average was no greater than 180, the gray value of the pixels on the *i*-th approximate ring was processed as follows, and the pixels in the whole image whose gray value exceeded 180 were made equal to 180:

$$
gamma\_{x,y} = \text{gra}\_{x,y} + 180 - \text{Average}\_i \tag{3}
$$

where *grayx,y* represents the pixels whose gray value was less than 180 on the *i*-th approximate ring; *x, y* indicate their horizontal and vertical coordinates; and 180 is the target gray value.

The one-dimensional entropy of an image reflects the average amount of information contained in the image, and provides the aggregation characteristics of gray distribution for the image. The larger the entropy is, the greater the probability of different levels of gray in the image. One-dimensional entropy of the image was introduced to evaluate the brightness homogenization process. The formula is:

$$H = -\sum\_{0}^{255} p\_i \times \log\_2 p\_i \tag{4}$$

where *H* represents the image entropy and *pi* indicates the probability of the *i*-th gray value appearing in the image, obtainable from the gray histogram.

The results of image entropy calculation before and after processing were between 5.336 and 4.869. The entropy value of images after brightness homogenization was shown to be significantly lower than before processing, suggesting that the gray value of images after processing was more uniform and the processing effect was satisfactory.

The problem of artifacts was effectively reduced after the image had been processed by brightness homogenization as shown in Figure 6. Additionally, the differences in grayscale between the central region and the edge region of the test piece were reduced significantly, and the brightness of the image was more uniform, which facilitated the threshold segmentation of the image.

**Figure 6.** Diagram of the effect of luminance homogenization. (**a**) Original image; (**b**) image after brightness homogenization; (**c**) grayscale contour map of original image; (**d**) after brightness homogenization.

#### 4.1.3. Image Filtering and Noise Reduction

The existence of noise makes the threshold segmentation program prone to errors in the recognition of aggregate edge, and can cause adverse effects on the calculation of threshold, as a result of which the threshold segmentation program shows sensitivity to noise. It was thus necessary to filter the image for the purpose of noise reduction. In this study, the adaptive low-pass Wiener filter [29] was applied to filter the grayscale images, with 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11 selected respectively as filter windows, to estimate the mean and standard deviation of local images, shown in Figure 7.

**Figure 7.** Comparison of filtering effects using filters of different sizes.

It can be seen in the figure that the amount of noise around the aggregate particles diminished with the increase of the filter window. Furthermore, the gray value of aggregate became more uniform. Meanwhile, however, the application of a larger filter caused severe loss of image details, thus resulting in deterioration of image clarity and aggregate edge identification. In this study, 7 × 7 filter window was applied to reduce the image noise.

The effect of the 7 × 7 filter is shown in Figure 8, revealing that the peak of the representative aggregate in the gray histogram narrowed but shifted upward after noise reduction, indicating that the aggregate gray value tended to be uniform and that the effect of noise reduction was evident.

**Figure 8.** Diagram of noise reduction effect. (**a**) Before noise reduction; (**b**) after noise reduction.

#### 4.1.4. Global Otsu Threshold Segmentation

Depending on the grayscale characteristics of the image, the Otsu algorithm divides the image into two parts: background and foreground, with the variance taken as a measurement of the uniformity of grayscale distribution. Thus, the greater the variance between the background and foreground, the more significant are the differences between the two parts of the image. In this situation, the threshold calculated by the Otsu algorithm is considered best when the variance between the two parts is the largest. Therefore, Otsu thresholding based on the whole image was selected to perform binarization. At this stage, the gray value of voids and artifacts in the image was returned to zero, and the image contained as few as three parts, i.e., black background, asphalt mortar, and aggregate. The aim of binarization is to separate the aggregate. Therefore, in this study, Otsu was applied

to calculate the double threshold. The smaller threshold was the boundary between the background and the asphalt mortar, while the larger threshold was the boundary between asphalt mortar and aggregate. The higher threshold was used for segmentation, with the gray value of mortar restored to zero. The result is shown in Figure 9.

**Figure 9.** Diagram of the adaptive processing method's threshold segmentation. (**a**) Before threshold segmentation; (**b**) after threshold segmentation.

#### 4.1.5. Image Morphological Processing

After thresholding segmentation, the holes require filling and conglutinated aggregates need to be separated. The holes distributed in the image can be classed into two categories. One results from the uneven texture of the aggregate and causes the low grayscale after CT scanning. This can be eliminated by the algorithm, and holes of this kind need to be filled. The other category relates to the proximity between the two adjacent aggregate edges, where the gaps between them become holes after threshold segmentation, and should be retained. This lays a foundation for the subsequent separation of the adhered aggregates by the watershed segmentation algorithm.

In this study, the image was first processed in reverse color, and then eroded [22], as a result of which the adjacent aggregate was separated. In the meantime, the holes inside the aggregates expanded, while remaining closed holes. The closed holes were filled, the holes inside the aggregates were removed, and the image was dilated to restore the size of the aggregates. Finally, the watershed algorithm [17,18,30] was applied to separate further the adhered aggregates. The outcome of the treatment is shown in Figure 10.

**Figure 10.** Comparison of before and after morphological treatment. (**a**) Before morphological processing; (**b**) after morphological processing.

#### *4.2. Ring Division Method*

With regard to the CT image processing of asphalt mixture, it has been proposed by some scholars that the mixture can be divided into multiple ring blocks, and that the impact of artifacts on threshold segmentation can be mitigated by threshold segmentation being performed in each part of the image. The final image can be obtained by combining the images after binary processing. The core idea is that the artifacts are distributed in a non-linear way along the radial direction of the test piece. In general, if segmentation based on the global threshold is adopted directly, the closer the gray value of the aggregate to the center of the test piece, the lower is the gray value difference. A more significant gray value difference in the aggregate will contribute to the removal of the aggregate from the center of the test piece. If the image is divided into ring blocks, however, the difference in the aggregate's gray level in each part of the image can be reduced, while the precision of threshold segmentation can be improved.

The asphalt-mixture CT image was split into three rings and one circular surface, in line with the ring division rule, with each part overlapping each other. Each part of the image was filtered and denoised, the threshold was segmented, and the combined image was morphologically processed, as shown in Figures 11 and 12, respectively.

**Figure 11.** Partition of original image.

**Figure 12.** Diagram of the image-processing method based on ring division.

It was found that the CT image-processing method based on ring division was effective for mitigating the impact of artifacts on threshold segmentation, by splitting rings. The loss of aggregate recognition was less severe than with global threshold segmentation, and boundaries were clearly seen when the binary image was combined with a complete image. Compared with the adaptive CT image-processing method, an evident gap was observed in the completeness of the aggregate recognition.

#### **5. Grading Verification of 3D Reconstruction Model**

To validate the image-processing method, the images processed using an adaptive image-processing method, the ring division method, and a global threshold segmentation method, respectively, were imported into Avizo software for 3D reconstruction [31,32]. The model is shown in Figure 13.

**Figure 13.** Three-dimensional reconstruction model of asphalt mixture.

Applying the extracted volume information for the particles in the model, the particles were regarded as spheres, and virtual sieving was conducted. For the virtual sieving, the software automatically identified the volume of the particles and equate each particle with a sphere of the same volume. The diameter of the sphere was considered to be the diameter of the particle. Because the quality of the CT images was not as high as required, the accuracy of aggregate extraction below 4.75 mm was far from satisfactory. Furthermore, it was difficult for some of the fine aggregate to be separated by the algorithm, due to its close contact and adhesion with coarse aggregates. Its small area caused some fine aggregates to be eroded. Consequently, when the grading curve was drawn, it was suitable only for counting aggregates with a diameter of 4.75 mm or over. The true grading required conversion to a diameter of 4.75 mm or more. The comparison of gradation is shown in Figure 14. When model gradation was evaluated, it was also possible to regard the model gradation as a regression fitting to the true gradation, and to use statistical parameters to evaluate the regression fitting effect.

As shown in Figure 14, it was found that the three-dimensional reconstruction model established from the images processed with the adaptive method was the closest to the real gradation. The error analysis is shown in Table 4. The error of the model using the adaptive method was smaller compared with the method based on ring division, and significantly smaller compared with the method based on global threshold segmentation, suggesting that the adaptive method is more accurate than other methods for the identification and segmentation of aggregates. Thus, it can restore the true gradation of the aggregates to the greatest extent. It was also revealed that there was little difference in accuracy of aggregate recognition between different kinds of asphalt mixtures, which suggests that asphalt mixtures of different kinds have good adaptability to the application of CT images.

**Table 4.** Gradation error analysis.

**Figure 14.** Comparison of three test pieces' gradation curves. Adaptive: image processing based on the adaptive image-processing method. Ring: image processing based on ring division. Global: image processing based on global segmentation. (**a**) SMA grading comparison; (**b**) OGFC grading comparison; (**c**) AC grading comparison.

#### **6. Conclusions**

curves and the real gradation curves were characterized by statistical parameters. It was discovered that compared with other methods, model gradation caused the least significant error from the real gradation, the degree of aggregate adhesion was also lower, and the adaptive image-processing method showed a strong adaptability to the processing of CT images for different kinds of mixtures.

**Author Contributions:** L.Z. (Lei Zhang) Writing—review & editing, G.Z., Y.T. Data curation, K.Z., C.X. Writing—original draft, Y.W. Investigation, C.C. Resources, L.Z. (Liting Zhao), X.L. Software, J.X. Validation, X.L. Visualization, L.W. Methodology, Y.T. Supervision, C.X. Funding Acquisition. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper was supported by the National Key Research and Development Program of China (No. 2018YFB1600203), National Natural Science Foundation of China (No. 51908168), Postdoctoral Science Foundation of China (No. 2019M651192 and No. 2020T130150), and Heilongjiang Postdoctoral Fund (No. LBH-Z19163).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the data available in main text.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **References**

## *Article* **Improved Procedure for the 3D Reconstruction of Asphalt Concrete Mesostructures Considering the Similarity of Aggregate Phase Geometry between Adjacent CT Slices**

**Chao Wang 1, Hui Xu 1, Yan Zhang 2, Yiren Sun 1,\*, Weiying Wang 3,\* and Jingyun Chen 1**

**\*** Correspondence<EMAIL_ADDRESS>(Y.S<EMAIL_ADDRESS>(W.W.)

**Abstract:** Existing image segmentation algorithms used for the computed tomography (CT) images of asphalt concrete mostly ignore the similarity of aggregate phase geometry between adjacent CT slices, thus increasing the variability in the aggregate phase pixel values between adjacent slices and leading to a large number of model defects, e.g., interconnected aggregates, flaky aggregates, and incomplete aggregates. The developed mesostructural models with these defects pose a challenge to following simulation operations. To address this issue, an improved procedure for the 3D reconstruction of asphalt concrete mesostructures considering the similarity of aggregate phase geometry between adjacent slices was developed, which includes two adjacent-slice pixel-value-correction algorithms, a multidirectional multiple-correction method, and an image pixel interpolation process. First, the bilinear interpolation algorithm was employed to improve the pixel density of 2D CT images and the average filtering algorithm was used to reduce the noise of the CT images. Subsequently, the OTSU method was employed to separate the asphalt mortar matrix phase from the aggregate phase, and the marker-based watershed segmentation method was used to separate the interconnected aggregates. Finally, the adjacent-slice pixel-value-correction algorithm was used to recover the similarity of aggregate phase geometry between adjacent CT slices, and the multi-directional multiple-correction method was used to further enhance the geometric similarity. The results show that the developed 3D reconstruction procedure removes most of the model defects in the 3D mesostructural model of asphalt concrete, thus realistically maintaining the 3D spatial distribution features and contour characteristics.

**Keywords:** asphalt concrete; 3D reconstruction; digital image processing (DIP); computed tomography (CT); mesostructure

#### **1. Introduction**

In traditional pavement analysis and design, asphalt concrete is often considered a homogeneous material and its mechanical properties can be determined using relevant laboratory tests [1,2]. However, asphalt concrete is essentially a heterogeneous material, which actually contains the aggregate phase, asphalt mortar matrix phase, and air void phase at the mesoscale [3,4], and thus its overall mechanical behavior is closely associated with the interaction between various components in the mesostructure.

Numerical simulation approaches, such as discrete element (DE) and finite element (FE) methods, have been applied to investigate the mesomechanical behavior of asphalt concrete; however, a reliable mesostructural model of asphalt concrete needs to be realistically established before using these numerical methods to deeply explain the mechanical behavior and damage evolution mechanism of asphalt concrete at the mesoscale. Currently, the methods for establishing asphalt concrete mesostructures are classified into two major categories, namely the random aggregate generation method [5–8] and the CT-based digital image processing (DIP) technique [9–11]. Compared with the random aggregate generation method,

**Citation:** Wang, C.; Xu, H.; Zhang, Y.; Sun, Y.; Wang, W.; Chen, J. Improved Procedure for the 3D Reconstruction of Asphalt Concrete Mesostructures Considering the Similarity of Aggregate Phase Geometry between Adjacent CT Slices. *Materials* **2023**, *16*, 234. https://doi.org/10.3390/ ma16010234

Academic Editor: Karim Benzarti

Received: 21 November 2022 Revised: 22 December 2022 Accepted: 22 December 2022 Published: 27 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the CT-based DIP technology can better reflect the real mesostructural morphological characteristics and spatial distribution of aggregates. Two types of CT-based DIP approaches, i.e., the 2D and 3D methods, are involved. The 2D methods are now more widely used in the asphalt paving community due to their simplicity in simulation. However, the 2D methods cannot reasonably explain the mechanical behavior because aggregates in asphalt concrete actually interlock with each other in a 3D fashion. In the common procedure for 3D CT image processing, filtering algorithms are first employed to reduce the noise of the 2D CT slice grayscale images. Then, image segmentation algorithms are adopted to segment the aggregate phase and the asphalt mortar matrix phase in 2D grayscale images to obtain 2D binary images. Finally, the processed 2D binary images of asphalt concrete are reconstructed into a 3D mesostructural voxel model along the CT acquisition direction. In this procedure, image segmentation is a key step that aims to facilitate the development of a geometrically reasonable 3D model by separating the aggregate phase from the asphalt mortar phase in 2D slices and eliminating interconnected aggregates in 2D slices. In the subsequent simulation, using these geometrically reasonable 3D numerical models will help improve the modeling accuracy and ensure computational convergence.

Among image segmentation techniques, the manual threshold selection (TH) [12] is the most basic one. The efficacy of this method is limited by the skill proficiency of the operator and it cannot guarantee the reproducibility of the image processing results. Therefore, it is not suitable for the high-volume image processing task required for 3D reconstruction [13]. OTSU is an automatic threshold determination method [14]. This method first assumes a threshold value in the grayscale range of image pixels, then classifies the pixels in the grayscale image into target and background categories according to the assumed threshold value, and finally uses the variance of the grayscale values of the two categories as the judgment index to obtain the best threshold value for this image by exhaustively searching for the maximum value of the variance [15,16]. Gong et al. [17] developed a 2D OTSU method. This method improves the segmentation quality of low signal-to-noise ratio images by simultaneously considering the gray values of image pixels and their 2D coordinate information [18]. Bhandari et al [19] developed the 3D-OTSU method. This method simultaneously considers the pixel gray value, neighborhood pixel gray median, and neighborhood pixel gray average to improve the segmentation quality of the image with a low signal-to-noise ratio and low contrast [20].

Although both the TH and OTSU methods can segment the aggregate phase and asphalt mortar matrix phase in the 2D CT slice images of asphalt concrete, a large number of model defects can be found in the obtained 3D model, such as interconnected aggregates, flaky aggregates, and incomplete aggregates, after 3D reconstruction is performed on the CT slice sets that have undergone the 2D image segmentation. Actually, there are huge differences in aggregate morphology for 2D slices and 3D models. The main reason for this discrepancy is that existing studies mostly ignore the connection of the pixels characterizing the aggregate phase between the adjacent slices in the procedure of obtaining optimal 2D binary image segmentation results, which inevitably increases the geometric differences in the aggregate phases between adjacent slices, and the accumulation of such geometric differences eventually leads to a large number of model defects in the constructed 3D voxel models. These model defects can adversely impact the convergence of simulation operations and even the correctness of simulation results. Due to the presence of numerous aggregates in the asphalt concrete mesostructural model, these defects are difficult to remove manually using software.

In the set of asphalt concrete CT slices, the 2D geometry and location information of the same aggregate possesses a natural similarity between adjacent slices. From a pixel perspective, this connection manifests itself as a continuity of pixel information at the same coordinates between adjacent slices. In order to overcome the above-mentioned model defects, this study developed a new procedure for the 3D reconstruction of asphalt concrete mesostructures considering the similarity of aggregate phase geometry between adjacent slices, which includes two adjacent-slice pixel-value-correction algorithms, a multi-directional multiple-correction method, and an image pixel density increase process.

#### **2. Methodology**

#### *2.1. Preparation of Asphalt Concrete Specimens*

In this study, asphalt concrete was considered as a two-phase particulate composite consisting of coarse aggregates (greater than 2.36 mm) and an asphalt mortar matrix containing the remaining fine aggregates, mineral fillers, and air voids.