#### Abstract

Rutting is one of the most important distresses for asphalt pavement structures especially due to the heavy traffic loadings. In this paper, the solution of Burgers model is employed and simplified, which makes the residual strain only depend on the applied load. Then, the governing equation for dynamic response of pavement was established and a semianalytical solution related to the stress state of asphalt pavement was obtained. Further, an explicit expression was presented and validated to calculate rutting depth of asphalt pavement through incorporating the dynamic stress solution, simplified Burgers model solution, and time hardening model. Compared with the FEM method, the calculation error is less than 5%, which proved the method proposed to be feasible. Compared with the numerical simulation method, every derivation of the rutting formula in this paper has the clear physical meaning and is hence easy to be understood. Besides, the proposed method can save a significant amount of time consumed due to its fast and reliable calculation.

#### 1. Introduction

Rutting is one of the most important distresses for asphalt pavement structures. It is caused by material consolidation and lateral movement due to repeated heavy wheel loadings on the various pavement structures. The distress is manifested by a depressed rut along the wheel path on the pavement surface. The rutting distress is viewed as not a structure failure but as a serious safety hazard to vehicles because hydroplaning can occur in the presence of rutting in rainy weather, resulting in serious traffic accidents [1, 2]. Moreover, vehicles tend to be pulled towards the rut path, making it difficult to drive. With the increase of road grade, construction of high-grade highways, improvement of traffic management, and emergence of channel traffic, rutting has become a critical issue for the design, construction, and use of asphalt pavement structure. Therefore, in addition to the high-temperature stability of asphalt mixtures, there currently is an urgent problem to be solved in the design, that is, how to predict the rutting of asphalt pavement structure from the mechanical properties of asphalt mixtures [3, 4].

Currently, a large number of indoor and outdoor studies on the rutting problem have been carried out by the international researchers. Many methods have been presented in the pursuit of a perfect design for asphalt pavement structure. In the earliest considerations of rutting in road design, the vertical strain on the top of the subgrade was limited to constrain the rutting [5, 6]. Those methods lacked a comprehensive approach and there was no direct analysis of the asphalt pavement layer, which has the largest effect on rutting depth. Some researchers performed creep experiments on asphalt pavement and studied the experimental data to analyze the rutting of asphalt pavement [7–13]. Although many indoor experiments had been conducted to simulate a real vehicle-road system under the experimental conditions such as environment and applied load, there are significant differences in the load system and the environment conditions between the indoor experiments and those encountered in practice. Due to the limitations of the experimental methods and simulation technologies, it is difficult to simulate the outdoor environment and the mechanical load effect. Hence, predictions of rutting from indoor experiments are still different from the actual deformation due to rutting. Some scholars studied the rutting of asphalt pavement through rutting experiments on asphalt pavement and analysis of the experimental data [14]. The experimental methods are relatively intuitive and provide a good analogy of the rutting formation process of asphalt pavement. High-precision experimental equipment is not required. The operation methods are not complicated and are easily used in engineering. However, there are a few disadvantages of this method. The basic mechanical parameters of the materials cannot be obtained, the method cannot fully reflect the actual road performance of a structure or material, and the experimental results cannot be used for theoretical calculations in pavement design. Software using the finite element method (FEM) was also used to predict the rutting of highway pavement [4, 15–18]. The results have been widely accepted. However, FEM software calculations are complicated and time-consuming as it is necessary to divide the domain into cells and the number of cells affects the results of the calculation. Ramsamooj et al. [19] studied the plastic strain of asphalt pavement and used it to calculate rutting depth. Lu and Sun [20] researched rutting damage of the semirigid base layer of asphalt pavement in China. They analyzed the formation mechanism of the rutting and the main factors of rutting. With a semiempirical and semitheoretical analysis method, a predictive model was built for the rutting and was found to be suitable for the semirigid base layer of asphalt pavement used in China. They considered the overall effect of the pavement structure in the semiempirical and semitheoretical analysis method. In comparison with the empirical method, their method reduced some limitations on the scope of use to a certain extent. However, their formula is still characterized by poor extrapolation. In addition, the application of mechanical theory in their work is not reasonable. Hence, there are still many difficulties in the application of this model.

It is well known that the Burgers model is well suited to simulate the viscoelastic mechanical behavior of an asphalt mixture [21]. The analytical solution of the creeping governing equation has been obtained. In the Burgers model, the strain is a function of the loading stress and loading time , while the recovery strain is a function of time. When the time tends to infinity, the recovery strain can be fully recovered accordingly. Then, the recovery strain depends only on the load stress. The residual strain is obtained by subtracting the recovery strain from the strain. When the load time is fixed, the residual strain depends only on the load stress. The residual strain can then be obtained by substituting the solved load stress into the simplified Burgers model. Hence, the premise is to obtain the residual deformation of the pavement by calculating the stress distribution of the load. The residual deformation of the pavement can then be obtained by using the analytical solution of the Burgers model.

In this work, depending on the characteristics of the semirigid base layer of asphalt pavement, a governing equation of the dynamic pavement response was established based on the theory of viscoelastic layered systems and the constitutive equation of soil. Based on the boundary conditions of the related hypotheses and the continuity of displacement and stress, the stress field distribution of asphalt pavement was obtained from the analytical solution of the governing equation [22]. Then, the solution of the Burgers model was simplified to calculate the residual strain. Meanwhile, a time hardening model for the asphalt pavement was introduced to define a predictive formula for rutting under multiple repeating loads. The accuracy of the model was verified by comparing the results from a numerical simulation with those from analysis software.

#### 2. Simplification of the Burgers Model Solution

Asphalt mixture is a kind of heat-flow material. In order to describe the thermal rheological properties of asphalt mixtures, Saboo and Kumar [23] applied the rheology theory for studying the dry and elastic characteristics of asphalt concrete. The Burgers model was proven to be suitable for analyzing the viscoelastic properties of asphalt mixtures [21] (Figure 1).

The deformation in the Burgers model includes three parts, namely, elastic deformation, viscous flow, and viscoelastic deformation (Figure 1).

The constitutive equation of the Burgers model is given by

Equation (1) can be rewritten as

Equation (2) is the constitutive equation of the Burgers model. It stands for a four-parameter fluid and can be used to describe the viscoelastic behavior of materials. This equation can be used to describe the initial two stages in a material’s creep curve.

The strain-time curve of Burgers model is shown in Figure 2 at , a constant stress is suddenly applied and maintained in place. The resulting strain can be obtained from the following:where .

If the load is removed at , the recovery strain iswhere .

For a highway, the time during which a car passes a certain location is very short. The interval between two cars passing by the same point is relatively larger than . Equation (5) can then be simplified by letting to completely recover the recovery strain:

The residual strain can then be determined by

It can be seen in (7) that parameters , , , , and are constants. The residual strain is a function of the load stress . Hence, once the load stress is obtained, the residual strain can be obtained as well. In the next section, the solution of the stress of the surface layer will be illustrated.

#### 3. Stress State of Asphalt Pavement under Moving Traffic Load

In this paper, the pavement system is treated as a two-dimensional, eight-layer system. The eight layers are the (upper, middle, and lower) surface layers, the (upper, middle, and lower) base layers, and the (upper and lower) subgrade. The surface layers, the base layers, and the subgrade are regarded as dry elastic media. The governing equation of the dynamic response of the pavement is established by ignoring the compressibility of the solid particles of asphalt mixture. The governing equations of the base layers and subgrade are represented by a dynamic equilibrium equation. The external load on the pavement surface is a moving bar-type load, which is used to express the moving traffic load [22]. The load expression is then converted into a Fourier calculation expression and substituted into the governing equation, which is then transformed from a partial differential equation to an ordinary differential equation for simplification. Based on the related hypotheses and the boundary conditions pertaining to the continuity of displacement and stress, the governing equation is solved with a semianalytical solution.

As shown in Figure 3, the pavement system consists of uniform elastic medium surface layers, base layers, and subgrade. The subgrade is fixed on a rigid base bed at a certain depth. Hence, the displacement at the rigid road bed is regarded as zero. Each structural layer of the pavement system is regarded as infinite laterally.

In order to establish the model of dynamic response in pavement, some assumptions are made as follows:(1)Each structural layer of the pavement is a medium of homogeneous, isotropic, and elastic material due to the transient load in a very short time.(2)The deformation of the asphalt concrete is very small.(3)The aggregates in asphalt concrete is incompressible.(4)The displacement and stress between the structural layers are continuous.

No consideration will be made of the gradual compaction process of the pavement and the shrinkage process of the pores. Without taking water into account, the following governing equations are obtained [24]:In (8), and are the positive stresses of the solids along the and directions, respectively (units: Pa), and is the shear stress along the plane (units: Pa).

For the two-dimensional planar strain problem, the following expressions are obtained from Hooke’s principle of stress and strain and the damping property of the material.

On the right side of (9), the minus sign designates that the value of a compressive stress is positive. and are the effective stresses along the horizontal () and vertical () directions, respectively; and are the strains along the and directions, respectively; and are the displacements along the and directions, respectively; is the Poisson ratio; and is the dynamic shear modulus of the solid material [25], where is the damping coefficient of the material, and is the static shear modulus of the solid material.

For the pavement layers, the base layers, and the subgrade, the dynamic governing equations are shown in (14), where is a layer indicator from 1 to 8, designating the upper surface layer, the middle surface layer, the lower surface layer, the upper base layer, the middle base layer, the lower base layer, the upper subgrade, and the lower subgrade, respectively.

##### 3.1. Expression of Moving Traffic Load

A moving load function () may be represented with a Fourier series as follows [26, 27]:In (16), the load has a width of and a strength of ; , where is the distribution period of the moving load. Here, ; and is a real number in the range from −∞ to +∞.

When using a moving coordinate system, that is, assuming that the coordinate system moves along with the moving direction of the load, then the expression of the Fourier series of the moving load of the traffic at any moment is

According to the Fourier transform, is

##### 3.2. General Solution of the Dynamic Governing Equation

For a linear system, an arbitrary function can be expressed as the sum of a series of linear harmonic functions [28], which can be written asIn (19), is a function of the single independent variable . It is an th harmonic function, while the exponential function is only a function of .

According to formula (19), the governing equation of the entire system can be written in the following form:The parameters in the above equations are () as follows:

The general solution of the governing equation may be expressed as follows:After simplification,where , and () are the roots of the 4th-order polynomials inThe solutions (roots) can be obtained through numerical methods (e.g., MATLAB program).

*Boundary Conditions and Solution*. In order to solve the governing equation, it is necessary to obtain the integration constant (.

At the top of the upper surface layer (), , .

At the bottom of the upper surface layer (), the boundary condition of the displacement is and . Then, the group of equations for the boundary condition can be written as

The boundary conditions for the middle surface layer, the lower surface layer, the upper base layer, the middle base layer, the lower base layer, and the upper subgrade () are

In the lower subgrade model, the boundary conditions are listed as follows.

At the interface between the lower subgrade and upper subgrade (), the boundary conditions for the displacements are , .

At the interface between the lower subgrade and upper subgrade (), the boundary conditions for the stress are , .

At the bottom of the lower subgrade (), the boundary conditions of the displacements are , . Then, a group of equations for the corresponding boundary conditions can be written as

For all the groups of equations of the boundary conditions shown above, the numerical methods can be used to obtain the displacement ( and ), the positive stress (), and the shear stress (). These solutions can then be substituted into the general expressions for the corresponding dynamic governing equations of the pavement system under a moving load. The solutions for the displacements can be expressed as the following form:

In addition, the expressions of the shear stress and the positive stress are

Therefore, the stress in pavement system under moving traffic load can be expressed according to the above derivation. Further, it is significant to validate the solution of stress state in pavement which can be illustrated as follows. Here, a FEM simulation is conducted, and the calculation parameters of pavement structure are listed in Table 1 [29].

Both the proposed method and FEM simulations are used to calculate the shear stress, vertical positive stress, and horizontal positive stress at the bottom of the asphalt pavement. The comparison results are shown in Figure 4.

It can be seen in Figure 4 that the calculation results (shear stress, vertical stress, and horizontal stress) of this work excellently match that of the FEM model. Therefore, the stress obtained through the presented solution can be used to evaluate the deformation in the following analysis.

#### 4. Establishment of Rutting Calculation Model

##### 4.1. Rutting Deformation

From the obtained stress field distribution of the pavement, it can be seen that the vertical stress of the surface layer is a function that fluctuant periodically. The action distance (stress at the action point begins to appear within the range of , as shown in Figure 3) is [m] at both sides of the observation point. The action time of the stress is , where is the moving speed of a vehicle. During the actual action time which is very short, the stress is fluctuant periodically as well; nevertheless, in the solution of Burgers model, the action stress is kept constantly during ( is load time, as shown in Figure 2). Therefore, needed to be transformed into during action time . According to the stress equivalence principle, can be given by

Substitute (32) into (7) and integrate the residual strain along the depth direction (the thickness of the surface layer) to obtain the following rutting depth for one load cycle, that is,

##### 4.2. Time Hardening Model

As the axle load number and time increase, the asphalt pavement will be gradually aged, and the viscous coefficient of asphalt pavement will gradually increase [30]. That is to say, the rutting depth cannot develop at the same speed of one load cycle (33). Thus, a time hardening model is chosen for restricting the rutting increasing speed; that is, the coefficient of viscosity of asphalt pavement is given bywhere is the number of standard axle load cycles and is related to the parameters of the asphalt pavement. Then, the rutting expression under a number of load cycles can be given by

In (35),

#### 5. Model Validation and Parameter Analysis

##### 5.1. Model Validation

Results of experiments on asphalt pavement described in An’s research [31] are used herein as an example to verify our analytical solution, and the FEM with ABAQUS [29] is also adopted to verify the calculation results. The parameters of elasticity and viscosity are evaluated and tabulated in Table 2 [31]. The operating temperature is 60°C and the corresponding deformations measured from the test, analytical method, and FEM are shown in Figure 5. Good matches can be seen between the analytical and FEM prediction and test data under the given ALN levels.

Furthermore, to extend comparison under various temperature, applied load, large axle load number, and so forth, the theoretical calculations of rutting for two kinds of pavement structures (Figures 6 and 7) are compared with the results of FEM modeling. The effects on rutting are comparatively analyzed for the factors such as the axle load number (ALN), axle load (wheel pressure), pavement temperature, and vehicle speed of vehicle.

The parameters of the asphalt mixtures in the Burgers model are shown in Table 3 [32], and the parameters of the nonasphaltic layers in this work are shown in Table 4 [32]. The results are summarized in Table 7. For the calculations of the two kinds of asphalt pavements, the parameter of the time hardening model () is equivalent to −0.42; the coefficient of rolling friction is 0.02.

It can be seen in Table 5 and Figure 8 that the calculation results of this work agree with that of the FEM model. The slight errors may be come from the model assumption which claims that the recovery time of pavement deformation tends to be infinite. Actually, the residual strain of asphalt pavement in a finite time does not reach the final value calculated by (7). Nevertheless, the strain can almost recover at finite time. In addition, the simplification of stress equivalence based on (32) may also lead to the difference between FEM and analytical results. However, the errors of asphalt pavement rutting in the calculation example are almost below 5%, which validates the method proposed in this paper. Furthermore, both the analytical and FEM methods are used to calculate the viscous flows in the rutting deformation of the asphalt pavement in the following analysis.

**(a)**

**(b)**

##### 5.2. Effect of ALN

Figure 8(a) shows the relationship between the RD and ALN for two pavement structures. Obviously, the power function appears. That is to say, RD increases with ALN in a power manner. It is worth noting that the rutting depth increased rapidly with the increase of the ALN at the initial stage, and then the increase of rutting decreases sharply and gradually keeps stable at a small value. Hence, the hardening rate of pavement material is relatively larger at the initial stage of traffic on the asphalt pavement and rapidly decreases with ALN (Figure 8(b)). For instance, when the ALNs are and , the calculated value of RD is approximately 3.1 mm and 6.4 mm, respectively. Nearly 50% RD can be reached under the condition of one-sixth ALN. As far as the application is concerned, the relation between RD and ALN can effectively predict the trend of the rutting depth and provide the information that more attention should be paid to the early period after the pavement is opened to traffic in the preventive maintenance of the highway.

##### 5.3. Effect of Wheel Pressure

The effect of wheel pressure on the RD is investigated, and wheel pressure varies from 0.7 MPa to 1.2 MPa, which are commonly applied on the asphalt pavement. The temperature is 60°C, vehicle speed is m/s, the coefficient of rolling friction is , and ALN is 1,000,000. The proposed method in this work and ABAQUS are employed to calculate RD of asphalt pavement structure A with the results shown in Table 6 and Figure 9.

From Table 6 and Figure 9, it can be seen that the RDs predicted by these two methods match well, and the maximum error is less than 5%. The RD has an approximately linear relationship with the wheel pressure. Significantly, The RD value seems sensitive to the wheel pressure. For instance, the value of RD linearly changes from 8.52 mm to 23.04 mm when the wheel pressure ranges from 0.7 MPa to 1.2 MPa. That is to say, the magnitude of RD will grow up with 170% when wheel pressure increases with 71%. It should be point out that this kind of increasing relationship is only based on the rutting analysis for the given parameters. It does not mean that the RD will increase to an unlimited value for a large wheel pressure because the asphalt pavement will fail. Furthermore, the result indicates that the heavy traffic load (wheel pressure) will expedite formation of asphalt pavement rutting. In the design stage of the highway pavement, the effect of overload on the pavement surface should be considered in terms of heavy haul vehicles.

##### 5.4. Effect of Pavement Temperature

As we know, temperature is a very significant factor which plays a role in affecting rutting deformation. It is mainly because the properties of asphalt pavement partly, especially the viscoelasticity parameters, depend on the environmental and pavement temperature (Asi, 2007; Motamed and Bahia, 2011). For example, the elastic modulus and viscosity coefficient of asphalt pavement will decrease with the temperature increasing. Therefore, the effect of pavement temperature on RD is conducted. For vehicle speed of m/s, wheel pressure of Pa, the coefficient of rolling friction , and , the calculation results of RD for different temperatures of 15°C, 40°C, and 60°C are shown in Table 7.

In Table 7 and Figure 10, the RD increased with the growth of temperature. The RD appeared to be linearly increasing from 2.65 mm to 4.92 mm when the pavement temperature increased from 15°C to 50°C. However, the growth of the RD is soaring at an increasing speed which turned to be a nonlinear form when the pavement temperature is greater than 50°C. In high-temperature area, the temperature of asphalt pavement often changes from 15°C to 60°C. The line slope between 50°C and 60°C is sharper than that between 15°C and 50°C in Figure 10. That is to say, the effect of high temperature on RD is evident compared with that of low temperature. Therefore, the special attention should be paid to the material design to prevent formation of rutting under high-temperature condition.

##### 5.5. Effect of Vehicle Speed

For temperature of 60°C, wheel pressure of Pa, the coefficient of rolling friction , and , the calculation results for the rutting depth of asphalt pavement at different vehicle speeds are shown in Figure 11.

Figure 11 shows the relationship between vehicle speed and RD. It can be seen that RD decreases nonlinearly with vehicle speed increasing overall. Figure 12 shows the percentage of decline in rutting depth at different speeds range. It can be seen that this kind of decreasing trend is gentle with the vehicle speed increasing. Conversely, the rutting depth will be more evident when the vehicle speed of vehicle becomes slow. In practice, at some special road sections such as road section with large longitudinal slope and intersection of roads, the vehicle travels at a relatively low speed which makes loading time relatively increased. As we know, the material of asphalt pavement is viscoelastic and depends partly on load time. Therefore, the RD develops accordingly.

#### 6. Conclusions

The solution of Burgers model is employed and simplified, which makes the residual strain only depend on the applied stress. Then, the governing equation for dynamic response of pavement was established and a semianalytical solution related to the stress state of asphalt pavement was obtained. An explicit expression was presented and validated to calculate rutting depth of asphalt pavement through incorporating the dynamic stress solution, simplified Burgers model solution, and time hardening model. Compared with the FEM method, the calculation error was less than 5%, which proved the method proposed to be feasible. Through parameters analysis, the following additional conclusions can be drawn.

The axle load number (ALN) makes the rutting depth increase in a power function manner and the rutting depth (RD) increased rapidly with the increase of ALN at the initial stage. As far as the application is concerned, the relation between RD and ALN can effectively predict the trend of the rutting depth and provide the information that more attention should be paid to the early period after the pavement is opened to traffic in the preventive maintenance of the highway.

RD is sensitive to the magnitude of wheel pressure and increases approximately linearly. It indicates the heavy traffic load (wheel pressure) will expedite formation of asphalt pavement rutting. In the design stage of the highway pavement, the effect of overload on the pavement surface should be considered in terms of heavy haul vehicles.

The RD appeared to be linearly increasing when the temperature is less than 50°C. However, the growth of the RD is soaring at an increasing speed which turned to be a nonlinear form when the pavement temperature is greater than 50°C. The effect of high temperature on RD is evident compared with that of low temperature and the RD of asphalt pavement will be more evident for slow speed of vehicle.

The solution presented in this paper can provide theoretical technical support for predicting rutting of asphalt pavement and can be beneficial for dealing with the rutting of asphalt pavement. In practice, the attention should be paid to the factors such as wheel pressure of vehicle, ALN, pavement temperature, and vehicle speed to prevent and delay the rutting formation.

#### Notations

: | Loading time of the stress (s) |

: | Stiffness parameter of spring 1 of the Burgers model (Pa) |

: | Stiffness parameter of spring 2 of the Burgers model (Pa) |

: | Coefficient of consolidation (-) |

: | Interstitial hydraulic pressure (Pa) |

: | Viscosity parameter of damper 1 of the Burgers model (Pas) |

: | Viscosity parameter of damper 2 of the Burgers model (Pas) |

: | Strain of spring 1 of the Burgers model (-) |

: | Strain of spring 2 and dashpot 2 of the Burgers model (-) |

: | Strain of dashpot 1 of the Burgers model (-) |

: | Strain in the plane (-) |

: | Stress (Pa) |

: | Poisson’s ratio (-) |

: | First derivative of stress (-) |

: | Second derivative of stress (-) |

: | First derivative of strain (-) |

: | Second derivative of strain (-) |

: | Elasticity modulus (Pa) |

: | Volumetric strain (-) |

: | The recovery strain (-) |

: | The residual strain (-) |

: | Moving traffic load (Pa) |

: | The distribution cycle of moving traffic load (m) |

: | The velocity of moving traffic load (m/s) |

: | Moving time (s) |

: | The position of moving traffic load (m) |

: | The width of moving traffic load (m) |

: | Imaginary component (-) |

, : | Velocities of solid particle in and directions, respectively (m·s^{−1}) |

: | Angular velocity (s^{−1}) |

: | Mass of solid particles (kg) |

: | Shear modulus (Pa) |

: | The dynamic shear modulus of the solid material (Pa) |

: | Poisson’s ratio (-) |

, : | Total normal stress in and directions, respectively (Pa) |

, : | Effective normal stress in and directions, respectively (Pa) |

, | The normal strains along and coordinates, respectively (-) |

: | Volume strain (-) |

: | Elastic shear stress (Pa) |

, : | Space coordinates in and directions, respectively (m) |

, : | The displacements along the and directions, respectively (-) |

: | Damping ratio of material (-) |

: | The dynamic shear modulus of the th pavement layer (Pa) |

: | Poisson’s ratio of the th pavement layer (-) |

, : | The displacements along the and directions of the -th pavement layer, respectively (-) |

: | Integer number (-) |

: | Distance of the th pavement layer bottom to top of surface course (m) |

: | Elastic shear stress of the th pavement layer (Pa) |

, : | Total normal stress of the th pavement layer in and directions, respectively (Pa) |

: | The roots of the 4th-order polynomials in (24) (-) |

The number of standard axle load cycles (-) | |

: | Related to the parameters of the asphalt pavement (-). |

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research has been supported by the National Natural Science Foundation of China (Grants nos. 51248006 and 51308554), the Special Financial Grant from the China Postdoctoral Science Foundation (Grant no. 2013T60865), and the Guizhou Transportation Science and Technology Foundation (Grant no. 2013-121-013) to the corresponding author. The research is also supported by the Guizhou Science and Technology Department Foundation (Grant no. 20132167), the Hunan Transportation Science and Technology Foundation (Grant no. 201237), and the Central South University Graduate Student Innovation Fund Project (2015zzts062).