Global existence of positive and negative solutions for IFDEs via Lyapunov-Razumikhin method

1.   Introduction

With the rapid development in the past two decades, China's total highway mileage has reached number one in the world. At present, a large number of asphalt pavements have reached the end of their design life. The maintenance, repair and renovation of old pavements have become essential to pavement design and construction projects. Cold recycled technology can make full use of the old material milled from old pavements and achieve the mixing of asphalt and aggregate under normal temperature conditions, thus reducing the consumption of aggregate resources and energy ^[1].

Research on cold recycled mixtures has been relatively abundant, covering the performance design of emulsified asphalt and foam asphalt ^[2,3,4], mix grading design ^[5,6], the effect of old material admixture on mix performance ^[7], early strength ^[8,9], high and low-temperature stability ^[10,11,12], water stability ^[13] and fatigue resistance of cold recycled mixtures ^{[13,14,15,16,17,18]}. The intermingling effect of emulsified asphalt and cement determines the early strength of cold recycled asphalt mixtures. The hydration products of cement are distributed in microstructure in the form of needles and clusters, which are interwoven with the asphalt film to form a spatial three-dimensional mesh structure wrapped around the aggregate. Still, their internal microscopic defects reduce the bonding force between asphalt and cement and the strength of the hydration products ^[19,20]. On the other hand, the agglomeration of the old material particle clusters and the complex interfaces within the mixture significantly influence the strength and deformation properties of cold recycled asphalt mixtures ^[12,21,22]. Studies on the high-temperature stability of cold recycled mixtures have shown that the dynamic modulus and phase angle of cold recycled mixtures with temperature and frequency show similar characteristics to those of hot mix asphalt mixtures, which still exhibit a time-temperature dependent stress-strain relationship at low stresses. The presence of cement hydration products can improve the dynamic stability of the mixtures and the dynamic modulus at high temperatures and low frequencies. It also reduces the viscous properties and improves the elastic properties ^[20,23,24]. In contrast, recycled asphalt pavement (RAP) content and asphalt aging degree significantly affect the high-temperature performance of cold recycled mixtures. Due to the impact of aging asphalt, with the increase of RAP, the elastic properties of recycled asphalt mixtures increase, and the viscous properties decrease, increasing dynamic modulus and dynamic stability ^[25]. The study on the fatigue performance of cold recycled mixtures showed that the fatigue performance of cold recycled asphalt mixtures is weaker than that of hot mixed asphalt mixtures. Its fatigue behavior is equally sensitive and dependent on its material composition ^[2,26,27].

Regarding material performance, cold recycled mixtures are generally weaker than hot mixed asphalt mixtures ^[8,14,15], limiting their application in pavement structures. The current Chinese specifications only specify the key technical indicators in the design, production and construction of cold recycled mixtures but do not recommend the structural form of cold recycled asphalt pavement. There is also a need for theories and principles to support the design of repair and maintenance structures. In addition, the current design method is mainly based on the linear elastic layer system theory ^[28]. The influence of temperature and other factors are considered primarily to introduce modified coefficients. This approach, summarized from the field data, can make up for the shortcomings of the mechanical theory, but it is very dependent on the data samples. For example, whether the modified coefficients fitted from the original structural forms still apply to a new structure type is still being determined.

This study aims to propose a simulation method that can couple the measured temperature and traffic volume and calculate the spatial and temporal distribution of fatigue damage. The viscoelastic model and traffic volume distribution coefficients were introduced to characterize the effect of temperature and loading conditions. Moreover, the damage variable was defined based on the fatigue equation to determine the structural fatigue damage. Considering the difference in the bearing capacity of the old pavement structure, the effect of the bearing capacity of base layers on the fatigue life of the cold recycled pavement structure will also be analyzed.

2.   Field experiments

2.1.   Original pavement structure condition

The test road on which the study is based is located in the Shandong Province section of the G22 expressway, and the original pavement structure and related conditions are shown in Figure 1. It can be seen that the original pavement structure has considerable and severe distress, and the pavement serviceability is poor. It can be observed from the drill core sampling at the distress that there are apparent penetration cracks inside the cement-stabilized aggregate (CSA) base layers. The bond between the CSA and the asphalt mixture layer has been destroyed, and the fracture contains noticeable sediment. The existence of reflective cracks provides a path for moisture to enter the interior of the pavement structure, resulting in the rapid damage development to the asphalt surface layer subject to vehicle loading, and finally presenting the form of damage coupled with multiple distress such as net-shaped cracking, rutting, sinkage and pumping.

When the top 15 cm surface layers were milled, the mechanical state of the old base layers and subgrade was tested using the methods such as falling weight deflector (FWD) and dynamic cone penetrometer (DCP). The back-calculated modulus of CSA bases and subgrade are shown in Figure 2. It can be seen that, compared with the initial modulus without damage, the damage of the CSA base layers in most sections is severe, and the average modulus is 1493 MPa, which is only about 20% of the design value. This indicates that most of the CSA bases have lost sufficient bearing capacity. In contrast, the bearing capacity of the subgrade is better, and the average value of its modulus can reach 66 MPa, which meets the requirements of the Specifications for the design of highway asphalt pavement (JTG D50-2017) for the subgrade modulus. It also verifies that the CSA bases can protect the subgrade well.

2.2.   The design and monitoring of the cold recycled pavement structure

2.2.1.   Cold recycled pavement structure and sensor burial scheme

Based on the above conditions, a cold recycled pavement test structure was proposed, as shown in Figure 3. The asphalt surface layers of the original design were milled. The distress existing in the CSA bases was repaired mainly by the following two measures: 1) for areas with significant net-shaped cracks, they were milled and replaced by new CSA materials; 2) for other cracks, the repair was carried out by penetrating cement mortar. Then, a 16 cm layer of emulsified asphalt cold recycled (EACR) mixture + 18 cm hot mix asphalt (HMA) surface layers, consisting of 4 cm SMA-13, 6 cm AC-20 and 8 cm ATB-25, were paved. For the EACR mixture, the optimum moisture content of the laboratory Marshall compaction is 4.7%; the amount of cement is 1.5%; the amount of emulsified asphalt is 3.2%; the additional new aggregate is limestone, and its content is 14%; and the rest is all RAP.

Considering the typical temperature-dependent characteristics of HMA and EACR mixtures, several temperature sensors were placed in the middle part of each layer to measure the variation of the temperature field inside the structure. Meanwhile, strain sensors were placed at the bottom of the ATB-25 layer and the bottom of the EACR layer to monitor mechanical responses located directly below the wheel track, two on each layer.

2.2.2.   Measured strain data

In collecting strain data, the test vehicle (it has a single rear axle similar to BZZ-100 but with a weight of 102.1 kN) passed over the section above the buried sensors at different speeds (5, 10, 25 km/h) to obtain strain response data for different layers; each speed will be repeated twice. The measured maximum strain value (tensile strain) and loading frequency of ATB and EACR layers and the structural temperature during the field tests are shown in Figure 4, where the load frequency was obtained by the inverse of the response time difference between the peak and the trough of the strain.

3.   Simulation model construction

3.1.   Model construction of pavement structure

Based on the cold recycled pavement structure, a three-dimensional pavement structure model was constructed in finite element software ABAQUS, as shown in Figure 5. According to the specification JTG D50-2017, a single-axle two-wheel standard axle load with an axle weight of 100 kN was used as the simulated cyclic load. It is generally simplified to a double-circular load of 0.7 MPa, 21.3 cm in diameter and 31.95 cm in the middle spacing. The loading speed is set to 60 km/h, which moved on the full size of the longitude direction (x-direction in Figure 5). Due to the setting of tack coat, seal coat and prime coat in field structure, the contacts between HMA-HMA, HMA-EACR, and EACR-CSA layers were considered fully continuous; and the contact between subgrade and CSA base layer is frictional, the friction coefficient is 0.3. To improve the simulation efficiency, the symmetry plane of the double-circle load is used as the reference plane to establish a symmetric model. The response was obtained from the middle part of the model. To improve the accuracy of the results, the mesh in this area and the loading track were densified. The model's dimensions are finally set to 10, 4 and 6 m in the longitude, transverse and depth directions by verifying the calculated results' stability. Symmetric boundary constraints are added to the symmetry surface of the model, and full constraints are given at the bottom to control the deformation and rotation in each direction.

3.2.   Traffic and temperature conditions

Considering the influence of calculation efficiency, the traffic axle load calculation only considers the monthly and hourly distribution in a day. First, the actual traffic volume of the test road for the whole year is counted, and all of them are equated to the equivalent cumulative cycles of the standard axle load according to the approach introduced in JTG D50-2017. Subsequently, the equivalent loading cycles for each month were counted, and the monthly distribution coefficient of traffic volume was obtained through Eq (1). At the same time, the hourly distribution coefficient of traffic volume can be obtained by averaging the hourly loading cycles for each day of the year, as shown in Eq (2).

where, ${M_i}$ represents the monthly distribution coefficient of the ith month; ${H_j}$ represents the hourly distribution coefficient of the jth hour; ${T_i}$ is the equivalent loading cycles of the ith month; ${T_{sum}}$ is the equivalent loading cycles of the whole year; and ${T_{k, j}}$ represents the equivalent loading cycles in the jth hour of the kth month.

The temperature field is constructed based on the measured temperature data. The day closest to the monthly average maximum and minimum temperatures is selected between the middle of each month (11th to 20th) as the representative day of the month. The temperature field under different hourly periods on that day is recorded as the characteristic temperature field of that month. Eventually, the temperature data recorded at different depths are linearly connected to obtain the continuous temperature field data. These can be written into the field variable definitions in ABAQUS through the UTEMP subroutine.

3.3.   Material constitutive models

HMA mixtures and EACR mixtures have typical viscoelastic properties, i.e., their modulus is correlated with the loading time and temperature. Their viscoelastic properties are characterized by the generalized Maxwell model as follows:

where $K\left( t \right), G\left( t \right)$ represent the bulk and shear relaxation modulus; and $\zeta \left( {\text{t}} \right)$ is the temperature adjustment function; it can be obtained by

where $a{}_T^{ve}$ is the temperature adjustment coefficient. In the thermal viscoelastic theory, the commonly used temperature adjustment coefficients are Arrhenius form and WLF form, where the WLF form ^[29,30] is selected in this study as shown in Eq (5):

where ${C_1}, {C_2}$ are material constants; and ${T_0}$ is the reference temperature.

The bulk and shear modulus can be expressed equivalently by Young's modulus $E$ and Poisson's ratio $\nu$ as:

For Young's modulus in viscoelastic models, the Prony series expression, which divides the dynamic modulus into moduli at several orders of relaxation time, is usually used to represent the generalized Maxwell model, as shown in Eq (7):

where $E\left( t \right)$ is the dynamic Young's modulus; ${E_e}$ is the static Young's modulus; ${E_i}$ is the modulus of the ith order; and ${\rho _i}$ is the relaxation time of the ith order.

The expression for the dynamic modulus in Eq (7) can be converted to an imaginary expression. Further, the parameters of each order of the Prony series can be obtained by fitting the stored modulus to the experimentally measured dynamic modulus and phase angle at different frequencies.

To directly obtain the fatigue performance of the cold recycled pavement structure, the damage variable D, defined by the ratio of current loading cycles to the maximum allowed cycles, is introduced to characterize the fatigue damage within the structure. The total fatigue damage for a specific service period can be calculated according to Miner's linear fatigue accumulation theory, as shown in Eq (8):

where ${d_{i, j}}$ represents the damage caused by the jth hour of the days in the ith month; and ${N_{i, j}}$ represents the maximum allowed cycles corresponding to the tensile stress ${\sigma _{i, j}}$ caused by the loading conditions of by the jth hour of the days in the ith month; and $a$ and $b$ are the material constants. The value of ${\sigma _{i, j}}$ uses the maximum tensile stress during one loading cycle; therefore, if the maximum value of the current calculated point is compressive, it will be considered as no damage accumulation occurred.

The above viscoelastic constitutive model and damage calculation process can be written into the UMAT subroutine of ABAQUS simultaneously, where the damage values are stored by the state variable STATEV so that the damage field distribution under repeated loading can be obtained.

4.   Numerical simulation

4.1.   Properties

The viscoelastic constitutive parameters of the HMA and EACR mixtures were determined by the measured dynamic modulus and phase angles at different temperatures and frequencies in the uniaxial dynamic compression tests. Eq (7) can be rewritten as:

where ${E_\infty }$ is the instantaneous Young's modulus. By introducing the modulus ratio ${e_i}$ in Eq (11), Prony's series can be expressed by Eq (12). Then the viscoelastic model parameters of the three HMA mixtures and EACR mixture can be determined as shown in Table 1.

For the cold recycled asphalt pavement structure, a certain level of damage has accumulated in the CSA bases during its service period due to the environmental effects of dry and wet, freeze-thaw cycles and traffic loading. Therefore, CSA bases' residual strength and stiffness are generally lower than newly constructed bases. However, results in Figure 2(a) show that the modulus dispersion of the CSA base in different road sections is large, which will inevitably affect the mechanical response of the whole reconstructed structure. To evaluate the influence of varying base states on the cold recycled asphalt pavement structure, the modulus of CSA bases is selected in the range of 1000–15,000 MPa. Meanwhile, the subgrade modulus is taken as the average of the test results 66 MPa, as shown in Figure 2(b). The constitutive model of both CSA and subgrade was a linear elastic model.

The measured traffic volume is processed according to Eqs (1) and (2), and the traffic volume distribution can be obtained as shown in Figure 6. To further simplify, by combining the traffic volumes for every two hours of a day, the traffic distribution under 12 time periods can also be obtained. According to the measured traffic volume, the annual total equivalent cumulative cycles of standard axle load for the test section is 2,648,440.

Similarly, the temperature data under 0:00, 2:00, 4:00, ... and 22:00 of each representative day of the month were selected and linearly interpolated to obtain the representative temperature field under these 12 time periods, which can be substituted into ABAQUS software for the calculation of mechanical response. Therefore, for each year of traffic loading, there are 12 × 12 loading groups, thus needing to divide into 144 analysis steps to calculate the mechanical responses and cumulated damage in turn. On this basis, the structure's service life can be obtained by introducing the annual growth rate of the traffic volume (which can be equivalently converted to the damage growth rate due to Miner's linear fatigue accumulation theory).

4.2.   Mechanical responses

According to the parameter values introduced in Section 4.1, the process of field tests can be simulated in the established finite element model. Table 2 compares the measured and simulated strain responses, where the loading speed was 10 km/h, and the temperature is shown in Figure 4(c). It can be observed from Table 2 that the strain response times and magnitudes of the bottom of the ATB layer and EACR layer are close to the data measured by the sensors, which also verifies the validity of the simulation model.

Figure 7 shows comparatively the stress response distribution of the pavement structure for old and newly constructed CSA bases under the temperature field at 12:00 in July. It can be seen that when the capacity of the CSA base is well, the stress distribution of the overall structure is closer to that of a conventional semi-rigid base asphalt pavement structure. The newly constructed CSA base, which means the modulus of the 43 cm CSA layers is equal to the newly constructed CSA, ensures that the bottom of the EACR layer is not in tension, making the CSA layer the main bottom-up crack-bearing layer. Due to the deeper position of the EACR layer compared with the ATB layer, the loading frequency of the EACR layer is lower, which, together with the material property difference, results in the lower modulus of the EACR layer than that of the ATB-25 layer, thus making a small amount of tensile stress appear at the bottom of the ATB layer. Also, it can be seen from Figure 7 that the location of the maximum tensile stress at the bottom of the ATB layer is close to the center of the loading circle, while when the layers are deeper (such as the EACR and CSA layers), the maximum tensile stress gradually transitions to the center of the double-circular load. For the old CSA base with severe damage, the modulus of which was set as 1500 MPa, the EACR layer's modulus is much higher than that of the CSA layer, leading to significant tensile stress at the bottom of the EACR layer. In this case, bottom-up cracks will most likely develop from the bottom of the EACR layer.

Figure 8 shows the stress distribution in winter. For the structure with the old CSA base, the modulus of both EACR and HMA mixtures will be higher than the modulus of the base layer due to the low temperature in winter and damage in the CSA base. It makes the bottom of the EACR layer subject to high tensile stress, which also becomes the layer where bottom-up cracks mainly develop. While in the structure with a newly constructed base, the high-strength CSA base can also protect the EACR layer and reduce tensile stress. However, since the low temperature makes the dynamic modulus of the HMA and EACR layers higher than that of the newly constructed CSA layer, the EACR layer can not appear to be free from tension under these conditions. The strain at the bottom of the ATB-25 layer fluctuated only in a low strain level range of 5.1 to 18.2 με throughout the year. From the study of fatigue properties of HMA mixtures and the design principle of permanent pavements, it can be assumed that this strain level will not lead to macro fatigue cracks.

Tables 3 and 4 demonstrate the strain and stress distribution of the EACR layer under different loading groups (time periods) at a base modulus of 1500 MPa. It can be found that the stress and strain responses of the EACR layer do not fluctuate significantly during one day due to the deeper position in which the EACR layer is located. On the contrary, the climatic temperature difference in different months is the main reason for the variation in the mechanical response of the EACR layer. Similar to the above analysis, the dynamic modulus of HMA and EACR mixtures become high due to low temperatures in winter, so there will be a large tensile stress at the bottom of the EACR layer. However, its strain level is relatively low owing to the high modulus. As reported in the study about the mechanical properties of EACR mixtures ^{[7,8,9,10,12,13,14,15,16,17,18,22]}, its splitting tensile strength is about 0.61–0.90 MPa. In this case, the results in Table 4 show that when the capacity of the CSA base is poor, the stress ratio at the bottom of the EACR layer can be close to about 0.5, at which the allowed fatigue cycles are very limited. Therefore, when the capacity of the CSA base is poor, it is not suitable to pave the EACR layer directly on the top of the CSA base layer; otherwise, in winter, it is very susceptible to fatigue cracking when subject to traffic loads.

4.3.   Damage calculation

In the structure shown in Figure 3, due to the thick asphalt layers (including HMA and EACR layers), the fatigue mode subject to traffic loads approaches the stress-controlled mode. Therefore, for the damage analysis of the EACR layer, the stress-controlled mode fatigue equation will also be used. Based on the previously described proportions, 100 mm diameter, 40 mm high cylindrical specimens were used for fatigue performance testing of EACR mixtures. After 28 days of curing at 25℃, the indirect tensile strength of EACR mixtures was first tested, and the average value obtained was 0.80 MPa. Then, the stress-controlled fatigue tests were conducted subject to the cyclic loads that were 40, 50, 60, 70 and 80% of the indirect strength. The cyclic loads were applied in the haversine form with a frequency of 10 Hz. The tests for each stress level were conducted four times, and the fatigue life of EACR mixtures are shown in Figure 9, which can be fitted by the Wöhler-curve as shown in Eq (13).

Then, taking the CSA base modulus of 1500 MPa as an example, the damage distribution in different areas of the EACR layer with different time periods can be calculated, as shown in Figure 10 and Table 5. It can be seen from Table 5 that the damage of the EACR layer shows the distribution characteristics of more in winter, less in summer, higher at noon and lower at night by the combined influence of temperature and traffic volume.

After obtaining the damage of the first year of traffic loading, the damage for each subsequent year can also be calculated by introducing the annual growth rate of the traffic volume based on Miner's theory. Further, fatigue life is the total service life of the pavement when the accumulated damage is equal to 1. The fatigue life of the EACR layer under different CSA base moduli were shown in Table 6. It can be observed that when the modulus of CSA bases gradually increased, the tensile stress at the bottom of the EACR layer gradually decreased, and its fatigue life also improved significantly (Table 6). However, in winter, the modulus values of both HMA and EACR mixtures will be considerably higher than those of the CSA due to the low temperature, leading to the inevitable damage accumulation at the bottom of the EACR layer. Therefore, no matter how to improve the base layer's modulus, the damage accumulation will still exist in winter, and the service life of the EACR layer is finite. When the modulus of CSA is between 9000–15,000 MPa (the modulus range of newly constructed CSA bases recommended by JTG D50-2017), it can be seen from Table 6 that the fatigue life of the EACR layer is about 6–14 years.

As the fatigue properties expressed in Eq (13), the fatigue resistance of the EACR mixture is weaker than that of the HMA mixtures, making it unsuitable to be used in the fatigue-bearing layer of the pavement structure. Thus, it needs to be controlled as much as possible so that the bottom of the EACR layer is not in tension during the structure design. Therefore, two main approaches are proposed to ensure the durability of the recycled structure. On the one hand, optimizing the fatigue resistance of EACR mixtures through material design is the most direct approach to improving the fatigue life of the structure. On the other hand, the structure durability can also be enhanced by transferring the fatigue cracking position from the EACR layer to other layers with better fatigue resistance. For example, a new CSA layer can be paved over the old base before paving the EACR layer. Or an HMA layer with higher asphalt content, which may have better fatigue resistance, can be added at the bottom of the EACR layer. In these cases, the tensile stress at the bottom of the EACR layer will be significantly reduced.

5.   Conclusions

This study proposes a finite element simulation method for cold recycled asphalt pavement structure that considers the effect of temperature variation and traffic volume distribution, which can quickly and effectively calculate the fatigue damage distribution of pavement structure and its service life. Based on the test road, the mechanical responses and fatigue damage of the cold recycled asphalt pavement structure with different CSA base modulus and different time periods were calculated by combining the field and laboratory test results. The main conclusions obtained are as follows:

1) The results collected by strain sensors show that the response frequencies of different pavement layers are mainly related to the depth where the layers are located.

2) Influenced by the temperature, the stress level at the bottom of the EACR layer shows a distribution pattern of high in winter and low in summer. Together with the influence of traffic distribution, the accumulation of fatigue damage in the EACR layer was more obvious in winter and at noon. Therefore, winter is the main period when fatigue damage will occur in the EACR layer for the test road structure.

3) The calculation results show that the relatively poor fatigue resistance of the EACR mixture will significantly affect the service life of the cold recycled asphalt pavement structure. Considering the large discrete modulus of old CSA bases, two main approaches were proposed to ensure structural durability: one is to optimize the fatigue resistance of the EACR mixture, and the other is to transfer the development location of bottom-up cracks.

A feasible framework for pavement damage calculation was proposed in this study, but there are still some limitations. For example, nonlinear damage evolution is not considered, and the damage behaviors of other materials and layers were not discussed in this study. Therefore, to improve the proposed method, the subsequent research will focus on the constitutive modeling of the pavement materials and field data observation and calibration.

Acknowledgment

The authors would like to acknowledge the financial support provided by Natural Science Foundation of Shandong Province (No. ZR2020KE024).

Conflict of interest

The authors declare there is no conflict of interest.